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

2008年汶川MS8.0地震前后地電阻率時空變化

2022-04-01 22:17:45王同利崔博聞王麗紅顏曉曄李菊珍
地震研究 2022年1期

王同利 崔博聞 王麗紅 顏曉曄 李菊珍

摘要:采用歸一化速率變化方法對2008年汶川MS8.0地震震中周邊19個定點臺站觀測到的地電阻率的時空變化過程進行了計算分析。時序分析結果顯示:在汶川地震前,震中周邊多個地電阻率臺站準同步地出現了中短期前兆異常;空間分布過程顯示,地震孕震過程中地電阻率在震中區域有負異常叢集和由遠及近的遷移現象,負異常叢集呈象限分布,長軸方向與震源機制解主壓力方向基本吻合,且受孕震斷裂、臺站位置、觀測裝置布設方向等的影響。

關鍵詞:汶川MS8.0地震;地電阻率;時序分析;空間分布;歸一化速率變化方法

中圖分類號:P315.722 文獻標識碼:A 文章編號:1000-0666(2022)01-0066-09doi:10.20015/j.cnki.ISSN1000-0666.2022.0008

0 引言

2008年5月12日14時28分,四川省汶川縣發生MS8.0強烈地震,造成69 000多人死亡,500多萬間房屋倒塌,有感面積達250萬km2,主震后又發生了多次余震。震后眾多的學者研究了可能與汶川地震孕震有關的前兆異常變化(杜學彬,2010;解滔等,2018;Luet al,2016;張學民等,2009;馮志生,2013;高曙德等,2010;何案華,2012;邱澤華等,2010)。

地電阻率是地下探測范圍內介質電阻率的綜合反映,其在地震孕震過程中表現出的異常變化主要反映構造應力作用下,測區介質變形誘發的微裂隙活動引起的電阻率變化,在靠近震中區的地電阻率異常形態常以負異常為主,并伴有年變畸變(杜學彬等,2007)。震前電阻率異常與震源機制有關,其升、降變化取決于臺址下介質的受力狀態(如膨脹或壓縮等),巖層電阻率對微小應變有放大效應,通??煞糯?個數量級。在地震發生前的短期階段,由于斷層預滑,近斷層破裂區的地下介質應變加速積累,地電阻率在原有的變化背景上通常會出現加速變化,而在近震中區應變被部分釋放,地電阻率通常表現為轉折的正異常變化(汪志亮等,2002;Du,2011)。對巖石進行應力加載實驗也表明,隨著應力的增大巖石電阻率呈下降變化,同時在相同應力加載作用下,巖石不同方向的電阻率呈現出各向異性變化(陸陽泉等,1990)。

我國的地電阻率觀測是在邢臺地震后快速發展起來的,通常采用對稱四極裝置,觀測供電極距在800~2 400 m,同一臺站在地表布設2~3個不同方向的測道,以整點觀測的方式進行長期連續觀測。在超過50年的觀測實踐中,多次記錄到了發生在臺網內和臺網附近的中強地震震中區的顯著地電阻率中短期異常(錢復業等,1982;錢家棟等,1985;桂燮泰等,1989;汪志亮等,2002;杜學彬等,2015;Lu et al,2016;解滔等,2018;王同利等,2020)。2008年汶川MS8.0地震前,其震中周邊大范圍內臺站記錄的地電阻率出現了較為顯著的異常變化(杜學彬,2010;Zhang et al,2009;Huang,2011;錢家棟等,2013;朱濤,2013;Lu et al,2016)。為研究汶川MS8.0地震前后地電阻率的時空動態變化,本文對其震中周邊19個定點臺站的地電阻率數據進行歸一化速率變化分析,研究地電阻率的時序變化過程和空間分布動態過程,探討地電阻率在強地震發生前后的異常變化、各向異性以及遷移過程等。

1 研究區概況

汶川MS8.0地震發生在青藏高原東緣的龍門山構造帶內,地震破裂面長達330 km。龍門山斷裂帶呈NE—SW向延伸,自新生代以來發生過多次強烈的地殼變形和破裂,致使該地區不同走向和規模的斷裂縱橫交錯,其西側是鮮水河斷裂帶,北側是西秦嶺北緣斷裂帶,西北側是東昆侖斷裂帶(圖1)。這些斷裂帶的地震活動頻度和強度均較高,其中鮮水河斷裂帶就曾多次發生7級以上大地震。

本文研究區范圍為四川及周邊的甘肅、云南、寧夏和陜西等地區(25°~40°N,95°~110°E),選取2003—2009年汶川MS8.0地震震中周邊500 km范圍內共19個臺站(圖1)記錄的地電阻率數據。臺間距為50~200 km,地表觀測極距為1 000~2 000 m,電極埋深為2 m;井下觀測極距為80~200 m,電極埋深為100~200 m。

2 研究方法

歸一化速率變化方法(Normalized Variation Rate Method,簡稱NVRM)主要用于處理原始時序數據曲線,消除年變化曲線等傳統方法不易識別的“弱”幅度異常(杜學彬等,2001)。計算過程為:首先對時序觀測數據中的個別數據進行預處理,然后再按固定步長ξ個數據組成子序列,對子序列進行線性回歸并求得回歸直線的斜率ki和相關系數Ri(第i個子序列),由ki·Ri得到初始速率值ρ·′i,以ξ為步長依次向后滑動計算初始速率時間序列;對初始速率進行均值為“0”“1”的歸一化處理,得到最終的時間序列(即ρ·值序列)。杜學彬等(2001)通過對196次MS3.2~7.9地震(94%以上為MS≥4.0地震)的計算表明前兆異常閾值為±2.4,無量綱。NVRM方法的關鍵是首先要對觀測資料進行預處理,去除突跳數據、年變化、非年變周期等,對存在階段性長趨勢變化的資料要根據實際異常落實情況進行分時段處理。該方法已得到廣泛應用,眾多學者獲得了地震前后的地電阻率時序異常(劉君等,2013;史紅軍等,2014;徐錫泉等,2014;閻睿等,2015;李姜,2019)。

NVRM原理如下,首先針對時間軸計算斜率:

ki=∑nj=1Tj∑nj=1yj-n∑nj=1Tjyj(∑nj=1Tj)2-n∑nj=1Tj2(1)

再計算自相關系數:

Ri=∑nj=1Tjyj-1n(∑nj=1Tj∑nj=1yj)∑nj=1Tj2-1n(∑nj=1Tj)2·∑nj=1y2j-1n(∑nj=1yj)2(2)

最后計算得出歸一化速率變化法序列為:

Si=Ri×Ki/σn-1 (i=n,n+1,…,N)(3)

式中:n為滑動步長;N為資料長度;σn-1為N-n個Ri×Ki時間序列的均方差,{y}是等間隔前兆數據時間序列,{T}是相應的等間隔時間序列。如果利用月均法曲線滑動,Si是月速率;如果利用五日均法曲線滑動,Si是五日變化速率。

3 地電阻率時序變化

對2003—2009年汶川地震震中周邊500 km范圍內臺站定點的地電阻率的原始觀測曲線進行分析,并進行NVRM分析,結果顯示:成都、冕寧、江油、天水、甘孜和武都臺的地電阻率長程時序曲線在震前均出現了中短期的下降變化異常(圖2)。圖2中橫虛線為歸一化閾值線,黑色區域表示超出閾值部分。

成都臺NE測道原始數據和NVRM計算數據曲線均在2006年中出現年變畸變,原始曲線出現持續性下降至地震發生后發生折返,下降幅度約5.5%(圖2a-1);NVRM曲線在2007年5月超過閾值,6月到達最低值后開始轉折回升,汶川地震就發生在折返上升的過程中(圖2a-2);成都臺NW與NE測道原始曲線變化基本同步,但變化較為不明顯(圖2b-1),NW測道NVRM曲線則更清晰地反映了震前地電阻率的變化過程(圖2b-2)。

冕寧臺EW測道原始曲線自2005年12月至2008年2月出現持續性下降變化,下降幅度約23%(圖2c-1),NVRM曲線顯示下降變化明顯超過閾值(圖2c-2)。NS和NW 測道2種曲線都無顯著變化。

江油臺NE測道原始曲線在震前無明顯變化,震后則出現快速下降-折返上升變化(圖2d-1),NW測道自2006年8月開始出現下降變化,下降幅度約1.1%(2e-1);兩測道的NVRM曲線顯示,在震前的下降變化均沒有超過閾值,在震后下降超過閾值后折返上升(圖2d-2、e-2)。

天水臺EW測道原始曲線于2008年2月底出現階躍性上升短期變化,上升幅度約1.5%,之后快速下降,汶川地震就發生在快速下降的過程中,地震后恢復至震前的觀測狀態,NVRM曲線變化基本同步(圖2f)。

自2006年開始,甘孜臺NE和NW測道原始曲線在震前2~3 a就開始出現持續趨勢下降,且年變形態清晰,震前累計下降幅度分別約為4.9%和8.3%,地震發生后原有變化趨勢也未改變(圖2g-1、h-1);NVRM曲線顯示,在2006年底,兩測道均出現下降變化,超過閾值后發生折返上升,于2007年5—6月上升到最高值后下降,同時NW測道年變形態發生改變,于2008年1月又開始出現下降變化,略折轉上升至3月再次下降,在汶川地震發生后又下降(圖2g-2、h-2)。

武都臺EW測道原始曲線年變形態清晰穩定,自2007年2月開始呈現趨勢上升,未出現年變畸變現象,地震后其變化趨勢也未改變,截至2009年6月上升幅度約8%,之后折返(圖2i-1);NVRM曲線顯示,震前上升超過閾值現象明顯(圖2i-2),地震發生在折返下降過程中,NS和NE測道震前無明顯異常變化。

4 空間異常動態分析

空間異常動態過程,是指異常的空間分布圖像隨時間遷移變化的過程。電阻率觀測的空間異常動態是由區域性應力場加強和介質不均勻性引起的,即除震源區是一個高應力集中區外,周圍介質還會形成多個集中點,這些應力集中點多分布在斷層附近,介質不均勻性較為明顯,且對震源區的應力變化較為靈敏。因此,當震源區應力場加強到一定程度時,周圍介質就會出現多個局部地段應力集中,形成區域應力場異常,這種反映震源區域的應力場變化通常會導致地電阻率的明顯異常,一般以地電阻率下降為主。強震發生時地電阻率下降的異常區范圍半徑約一二百千米。桂燮泰等(1987)、金安忠(1982)分析結果顯示在1976年唐山MS7.8地震前,華北地區相當大的范圍內地電阻率觀測值就出現了負異常隨時間、空間的動態變化??臻g上常常表現為以震中區為中心向外擴展,異常變化幅度從震中區向外逐漸由大變小,且分布不均勻;同一觀測站點的異常變化動態過程為隨著地震發生時間臨近,負異常變化幅度愈來愈大。王同利等(2020)分析了2017年九寨溝MS7.0地震前地電阻率的空間動態演化過程,結果顯示地電阻率負異常空間分布隨時間遷移,且和地震震中區有很好的關聯。

本文對汶川地震震中周邊19個臺站地電阻率觀測數據進行NVRM分析,得到了相對月變化,再進行空間等值線動態分析(圖3、4),計算中考慮大多數臺站的布極為NS和EW兩個測道,因此空間等值線動態分析為NS和EW兩個方向,布極方式沒有這兩個方向的臺站,分別用各臺和這兩個方向之間夾角最小的測道在NS和EW向上的投影計算得到的地電阻率變化量進行替代。

圖3顯示了2007年11月至2008年9月汶川地震震中周邊地區地電阻率NS測道的空間動態變化過程。從圖中可以看出,2007年11月震中區周邊地電阻率變化平穩;2008年1月在天水、通渭和武都附近出現了小區域的負異常變化;2月負異常變化范圍逐漸擴大,向NE-SW向延伸,SW向延伸至甘孜附近,NE向延伸至寶雞附近;3月負異常范圍進一步擴大且異常更加明顯,出現顯著的NE-SW向負異常條帶;4~5月負異常在震中區附近積聚,但相比3月有減弱的趨勢。5月12日汶川MS8.0地震發生,6月地電阻率負異常條帶現象減弱趨于消失,在7~8月出現近NS向的負異常條帶并向南遷移,2008年8月30日四川仁和發生MS6.1地震。綜上,隨著地震發生時間臨近,負異?,F象首先在震中區附近形成,然后面積逐漸變大并向外延伸形成條帶,且以震中區為中心呈象限分布,地震發生后異常條帶逐漸變弱消失或者發生遷移。

圖4顯示了2007年11月至2008年9月汶川MS8.0地震震中周邊地區EW測道地電阻率的空間動態變化過程。從圖中可以看出,2007年11月震中周邊地區地電阻率的空間變化基本平穩,沒有負異常現象;2008年1—4月在震中區附近出現由弱變強的NE-SW向負異常條帶,負異常區域首先出現在通渭和平涼附近,然后向NE和SW兩個方向延伸,NE向延伸至銀川,SW向經過成都、甘孜延伸至冕寧和紅格附近。2008年5月汶川地震發生后,負異常條帶減弱消失,6~9月僅有零散的負異常,大部分地區表現為正異?;蛘W兓?/p>

EW測道的地電阻率空間動態變化過程與NS測道基本一致,但異常變化稍顯明顯。EW測道負異常集中區域范圍大于NS測道,負異常強度也較NS測道強。地電阻率負異??臻g動態變化過程與地震孕育中的應力積累動態變化有關,且明顯地表現為各向異性。這種各向異性現象與九寨溝MS7.0、唐山MS7.8地震前地電阻率出現的負異常集中現象較為相似。汶川地震負異常集中區遷移方向、異常幅度增強方式與九寨溝MS7.0地震基本一致,這種異常遷移、震級不同、異常幅度也不同可能與地震斷層錯動方式和震中區巖性有關,兩次地震都發生在龍門山斷裂帶附近,汶川MS8.0地震的震源機制解顯示發震斷層以逆沖和走滑為主,震中區以花崗巖分布為主(易桂喜等,2012;王振榮,蘭江華,2008)。

5 結論

汶川MS8.0地震震中周邊地區地電阻率觀測是我國相對密集的地區之一。震中周邊共有19個地電阻率定點觀測臺站,多個臺站在地震發生前2~3 a出現了準同步下降變化。結合中強地震前地電阻率的孕震機理分析,獲得了以下主要結論:

(1)在汶川地震前,成都、冕寧、江油、天水、甘孜、武都臺地電阻率出現下降異常變化,從孕震機理上符合地電阻率地震異常變化機理,且臺站距離震中區域較近,這種變化是汶川地震前兆的可能性較大。

(2)地震前,震中周邊地區一個相當大的范圍內,地電阻率觀測值形成了以震中區為中心的電阻率負異常區域。異??臻g分布表現出比較顯著的條帶各向異性,長軸方向和震源機制解P軸、發震斷層破裂方向等基本一致。

(3)空間異常變化分布形狀取決于斷層方向、臺站和震中區的位置,并且與震中區附近的臺站數和臺站分布有關,異常變化幅度大小和臺站與震中區距離有關。

(4)汶川地震前地電阻率的負異常區域在NE—SW方向較為突出,負異常變化的長軸方向與震源機制解主壓力方向吻合,且EW測道條帶異常較NS測道明顯。這可能與在孕震晚期,近震中區主壓應力擠壓占優勢有關。地震前,孕震區域地下介質內部裂隙發育、裂隙走向沿主壓應力方向優勢排列,導電介質重新分布,導致地電阻率出現以快速下降為主的各向異性。不同觀測臺站或者同一觀測臺站不同測道的各向異性現象,是由于在孕震過程中震中區域應力加強引起發震斷層附近不均勻介質電阻率的變化所導致,與孕震斷層活動密切有關。

地電阻率時序分析能清晰地反映出汶川MS8.0地震孕震過程中斷層的應力變化過程??臻g動態演化過程分析是多臺觀測數據變化的集合,因而單臺數據受周邊環境的干擾影響可以忽略。時序分析和空間異常動態分析相結合,能夠更有效地揭示出中強地震的發生要素,對汶川地震的發震機理、震中區域的認識等有很好的幫助,但分析結果中難免存在因臺站密度低而造成偏差的情況。總之,汶川地震發生前,地電阻率出現長時程負異常變化、空間叢集現象真實地反映了強震孕育過程中發震斷層應力及周邊應力應變的變化,它和地震活動性等變化有很好的關聯性。分析地電阻率的時空變化過程對認識中強地震孕育過程及震源區應力場變化具有很大的意義。

參考文獻:

杜學彬,李寧,葉青,等.2007.強地震附近視電阻率各向異性變化的原因[J].地球物理學報,50(6):1802-1810.

杜學彬,劉君,崔騰發,等.2015.兩次近距離大震前成都臺視電阻率重現性、相似性和各向異性變化[J].地球物理學報,58(2):576-588.

杜學彬,阮愛國,范世宏,等.2001.強震近震中區地電阻率變化速率的各向異性[J].地震學報,23(3):289-297.

杜學彬.2010.在地震預報中的兩類視電阻率變化[J].中國科學:地球科學,40(10):1321-1330.

馮志生,李鴻宇,張秀霞,等.2013.地磁諧波振幅比異常與強地震[J].華南地震,33(3):9-15.

高曙德,湯吉,杜學彬,等.2010.汶川8.0級地震前后電磁場的變化特征[J].地球物理學報,53(3):512-525.

桂燮泰,關華平,戴經安.1989.唐山、松潘地震前視電阻率短臨異常圖像重現性[J].西北地震學報,11(4):71-75.

桂燮泰,關華平,戴經安,等.1987.地電視電阻率法測報強震指標討論[J].華南地震,7(2):56-63.

何案華,趙剛,劉成龍,等.2012.青海玉樹與德令哈地熱觀測井在汶川與玉樹地震前的異常特征[J].地球物理學報,55(4):1261-1268.

解滔,劉杰,盧軍,等.2018.2008年汶川MS8.0地震前定點觀測電磁異?;厮菪苑治鯷J].地球物理學報,61(5):1922-1937.

金安忠.1982.唐山地震前近震中區地電阻率的震前突變現象[J].地震學報,3(2):169-173.

李姜,喬子云,張國苓,等.2019.華北地區地電阻率歸一化速率異常分析[J].地震地磁觀測與研究,40(6):62-71.

劉君,杜學彬,范瑩瑩,等.2013.甘肅岷縣漳縣MS6.6地震前的地電阻率變化[J].地震工程學報,35(4):819-826.

陸陽泉,錢家棟,劉建毅.1990.大型花崗巖標本緩慢膨脹破裂過程中電阻率和聲發射前兆特征的實驗研究[J].西北地震學報,12(2):35-41.

錢復業,趙玉林,余謀明,等.1982.地震前地電阻率的異常變化[J].中國科學:化學,12(9):831-839.

錢家棟,陳有發,金安忠.1985.地電阻率法在地震預報中的應用[M].北京:地震出版社,10-51.

錢家棟,馬欽忠,李劭秾.2013.汶川MS8.0地震前成都臺NE測線地電阻率異常的進一步研究[J].地震學報,35(1):4-17,137.

邱澤華,張寶紅,池順良,等.2010.汶川地震前姑咱臺觀測的異常應變變化[J].中國科學:地球科學,40(8):1031-1039.

史紅軍,趙衛星,張可佳,等.2014.吉林省地電阻率震兆異常分析研究[J].地震學報,36(3):452-463,532.

汪志亮,鄭大林,余素榮.2002.地震地電阻率前兆異?,F象[M].北京:地震出版社,58-65.

王同利,崔博聞,葉青,等.2020.九寨溝MS7.0地震地電阻率變化時空演化分析[J].地球物理學報,63(6):2345-2356.

王振榮,蘭江華.2008.四川汶川大地震的構造分析[J].礦物巖石,28(2):1-10.

徐錫泉,高昌志,王亮.2014.內蒙古寶昌臺地電阻率長期觀測數據研究[J].地震工程學報,36(2):405-412.

閆睿,朱石軍,胡樂銀,等.2015.北京市延慶地震臺地電阻率長程觀測數據研究[J].地震工程學報,37(3):724-730,771.

易桂喜,龍鋒,張致偉.2012.汶川MS8.0地震余震震源機制時空分布特征[J].地球物理學報,55(4):1213-1227.

張學民,李美,關華平.2009.汶川8.0級地震前的地電阻率異常分析[J].地震,29(1):108-115.

朱濤.2013.汶川MS8.0地震前區域性地電阻率異常初步研究[J].地震學報,35(1):18-25.

Du X B.2011.Two types of changes in apparent resistivity in earthquake prediction[J].Science China(Earth Science),54(1):145-156.

Huang Q H.2011.Retrospective investigation oF geophysical data possible associated with the MS8.0 Wenchuan earthquake in Sichuan,China[J].J Asian Earth Sci,41(4-5):421-427.

Lu J,Xie T,Li M, et al.2016.Monitoring shallow resistivity changes prior to the 12 May 2008 M8.0 Wenchuan earthquake on the Longmen Shan tectonic zone,China[J].Tectonophysics,675:244-257.

Zhang X M,Ding J H,Sheng X H, et al.2009.Electromagnetic perturbations before Wenchuan M8.0 earthquake and stereo electromagnetic observation system[J].Chinese Journal of Radio Science,24(1):1-8.

Temporal and Spatial Evolution of the Apparent Resistivity beforeand after the 2008 Wenchuan MS8.0 Earthquake

WANG Tongli1,CUI Bowen1,WANG Lihong1,YAN Xiaoye2,LI Juzhen1

(1.Beijing Earthquake Agency,Beijing 100080,China)(2.Sichuan Earthquake Agency,Chengdu 610000,Sichuan,China)

Abstract

In this paper,the Normalized Variation Rate method is used to analyze the sequential and spatial variation of the apparent resistivity observed by 19 stations around the Wenchuan MS8.0 earthquake.The sequential analysis shows that some of the earthquake stations had quasi-synchronously recorded the short- and medium-term precursory anomalies of the apparent resistivity before the Wenchuan earthquake.The spatial analysis shows that in the epicentral region,the negative anomalies of the apparent resistivity had been featured with a strip distribution and with a migration from far to near before the earthquake.The cluster of the anomalies displayed a quadrant pattern,whose long axis was consistent with the direction of principal compressive stress of the focal mechanism solution.The seismogenic fault,the station location,and the layout direction of the observational instruments also influenced the distribution direction of the long axis of the precursory anomalies.

Keywords:the Wenchuan MS8.0 earthquake;apparent resistivity;time sequential analysis;spatial distribution;the Normalized Variation Rate method

主站蜘蛛池模板: 久久久久无码精品| 一级香蕉人体视频| 久热中文字幕在线| 亚洲欧美另类日本| 波多野结衣无码中文字幕在线观看一区二区 | aaa国产一级毛片| 国产成人综合亚洲网址| 国产精品大尺度尺度视频| 亚洲无码电影| 日韩免费毛片| 亚洲国产精品无码久久一线| 亚洲小视频网站| 国产亚洲欧美日韩在线一区二区三区 | 欧美中文字幕一区| 在线国产资源| 91啪在线| 91视频区| 毛片最新网址| 国产午夜不卡| 亚洲综合日韩精品| 国产三区二区| 亚洲大学生视频在线播放| 国产成本人片免费a∨短片| 一级毛片在线直接观看| 欧美一区二区自偷自拍视频| 91无码网站| 日韩精品免费在线视频| 久久精品国产国语对白| 免费一级毛片在线观看| 国产在线视频自拍| 一级片免费网站| 性网站在线观看| 国产成人精品综合| 91久久偷偷做嫩草影院电| 国产午夜一级毛片| 在线色国产| 亚洲三级网站| 美女无遮挡被啪啪到高潮免费| 日日摸夜夜爽无码| 992tv国产人成在线观看| 毛片免费在线视频| 亚洲精品日产精品乱码不卡| 在线日韩日本国产亚洲| 亚洲成a人片77777在线播放 | 无码免费视频| 九九精品在线观看| 日韩欧美视频第一区在线观看| 四虎精品黑人视频| 嫩草影院在线观看精品视频| 成人午夜天| 亚洲无码免费黄色网址| h视频在线观看网站| 无码福利视频| 曰AV在线无码| 丰满人妻久久中文字幕| a毛片免费观看| 香蕉久人久人青草青草| 丁香六月综合网| 国产福利免费观看| 欧美日韩中文字幕在线| 就去吻亚洲精品国产欧美| 亚洲无线国产观看| 日韩精品一区二区三区视频免费看| 91福利在线观看视频| 精品国产香蕉伊思人在线| 日本免费精品| 女人18一级毛片免费观看 | 亚洲精品动漫在线观看| 国产午夜不卡| 国产成熟女人性满足视频| 亚洲欧洲日产国码无码av喷潮| 欧美97色| 波多野结衣一区二区三区四区| 一级不卡毛片| 国产一二三区在线| 欧美国产日韩另类| 99久久精品免费观看国产| 无码人中文字幕| 日韩精品亚洲一区中文字幕| 亚洲天堂.com| 亚洲AV无码久久精品色欲| 久久77777|