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

川滇地塊現今地殼運動速度場解算與分析

2017-12-22 03:18:54徐克科伍吉倉
測繪工程 2017年3期
關鍵詞:特征

徐克科,伍吉倉

(1.河南理工大學 測繪與國土信息工程學院,河南 焦作 454000;2.同濟大學 測繪與地理信息學院,上海 200092)

?

川滇地塊現今地殼運動速度場解算與分析

徐克科1,伍吉倉2

(1.河南理工大學 測繪與國土信息工程學院,河南 焦作 454000;2.同濟大學 測繪與地理信息學院,上海 200092)

利用GAMIT10.4和BERNESE5.2軟件,分別對川滇地區中國地殼運動觀測網絡51個GNSS基準站2010—2014年的觀測數據進行處理。兩者解算所得單日解基線和坐標時間序列離散度優于10 mm,且變化趨勢一致,吻合度較高。基線重復率水平向優于1+4×10-9L;垂向優于4+6×10-9L(單位mm,其中L表示基線長度,單位m)。利用PL+VW噪聲模型進行了速度估計,獲取了川滇地區相對歐亞板塊地殼運動速度場。結果顯示:除滇西南地區呈現出由西向東的增大趨勢外,川滇塊體其它區域地殼運動特征由北往南、由西往東呈逐漸衰減趨勢;運動方向持續圍繞喜馬拉雅東結點作順時針旋轉,且仍然存在著南北向強烈的擠壓特征。以安寧河-則木河斷裂、小江斷裂和紅河斷裂為界,選用了兩側能反映震間斷裂構造特征的GPS測站,分析了斷裂兩側的速度差異。并利用平行斷層的速度剖面,擬合反正切函數變化趨勢,獲取了斷裂帶可能的變形寬度。

GPS;川滇地塊;地殼運動;速度場;斷裂帶

隨著中國地殼運動觀測網絡的開展,連續和分期GNSS觀測的持續,積累了時空分辨率越來越高的地表形變觀測數據。如何從一手的觀測資料中獲取高信噪比、真實可靠的測站位移序列或速度信息,為后續的地殼形變特征分析、斷層活動反演提供技術支撐。高精度的數據處理技術是重要的必備前提。通過對監測網絡基準站長年連續觀測得到的 GPS 數據進行精密處理和分析,可以獲得基準站的坐標時間序列,進而獲取研究區域高精度速度場。利用速度場數據通過區域動力學模型的正、反演數值模擬或計算,可進一步分析其地殼位移、應力、應變場等動力學特征,探討其構造運動的驅動力源和動力學演化模式。然而,同一軟件處理結果間的比較,僅說明了在一致的處理方案與原則下,計算結果的穩定性和內符合度。同種軟件處理結果間符合得再好, 不能完全斷定最終結果的無偏, 更不能估計系統性偏差的大小和誤差范圍。而不同軟件處理結果間比對,則可在一定程度上彌補以上缺憾,是一種比較客觀的評價方法。美國MIT的GAMIT/GLOBK軟件,瑞士伯爾尼科技大學的BERNESE軟件作為兩款最常用的GPS高精度數據處理軟件,這些軟件數據處理主要分為兩部分:一是對GPS原始數據進行解算獲得同步觀測網的基線解;二是對各同步網解進行整體平差和分析,獲得GPS網的整體解。盡管數據處理基本原理相同, 但在具體做法上和處理理念上差異明顯。兩種軟件結果的比較是檢驗網絡工程最終結果的一種有效手段。無論使用兩者中的哪一個,都可以利用另一個作為外部參考來進行驗證。只有這樣,才能客觀評價解算結果的可靠性。中國大陸川滇地區受青藏高原地殼物質的東向移動、阿薩姆頂點楔入及緬甸微板塊的共同影響,成為現今地殼運動最活躍的地區之一。而該地區活動斷裂的特殊構造特性, 又使得該地區的地殼運動表現為復雜多變的特征。為此,長期以來,川滇地區被視為重要的地震監測區域,并取得了大量的研究成果[1-6]。本文分別利用GAMIT/GLOBK10.4和BERNESE5.2兩款軟件對川滇地區2010—2014年間中國地殼運動觀測網絡基準站數據進行了解算,獲取了川滇地區相對歐亞板塊運動速度場,深入分析了現今川滇地區地殼運動特征和主要斷裂帶活動情況。

1 解算方案

為了能客觀地比較GAMIT/GLOBK10.4和BERNESE5.2兩種軟件的處理結果,盡量采用相同的處理參數,設置相同的參考框架ITRF2008、慣性系J2000.0、衛星截止高度角10°、IGS精密星歷和選用一樣的參考框架站等。

GAMIT/GLOBK10.4軟件解算方案:利用雙差觀測值,組成與觀測值和參數相關的非線性數學模型,整個GNSS網全部基線一起參與解算。解算中采用了最小二乘算法反復迭代來估計測站的相對位置、軌道、地球自轉參數、對流層天頂延遲參數和大氣水平梯度參數等,得到的載波相位整周模糊度分別為實數和整數的約束解與松弛解。選取了11個中國及周邊的IGS站聯合進行處理。本文參考框架點選擇的依據是:①全球均勻分布;②站點穩定,不受大震或大的構造活動的影響;③水平向重復性誤差不超過3 mm;④站點線性趨勢項明顯,一致性較好等。對IGS核心站設置強約束,水平向0.05 m,垂直向0.1 m?;鶞收静捎肂ERNESE PPP定位結果作為先驗值,水平向1 m,垂直向1 m。主要模型選擇與參數設置如下:衛星截止高度角設置10°,觀測值隨高度角定權。根據VMF1模型修正大氣層的映射函數。根據分段線性模型估計對流層天頂延遲和水平大氣梯度,每2 h估計1個天頂延遲參數。基線處理是利用偽距和載波相位觀測資料的雙差組合求得臺站坐標和衛星軌道的單日松弛解。由于基線向量無法提供確定點的絕對坐標所必需的絕對位置基準,GPS衛星的軌道和測站坐標都不在同一個穩定的參考框架里,整個GPS網作為一個剛體每次解算都會存在整體平移和旋轉。因此,必須引入外部基準。先利用雙差數據進行最小二乘參數估計,解算出各時段的基線和模糊度,然后將各同步觀測網自由基準的法方程矩陣進行疊加,再對平差系統給予確定的基準,獲得最終的平差結果。采用卡爾曼濾波的模型,對同步網解進行整體處理,獲取的測站坐標和速度。利用GLOBK將整個監測網每天的單日松弛解和圣地亞哥海洋研究所軌道中心SOPAC計算出的全球IGS跟蹤站的多個單日松弛解(igs1~igs6)合并,得到一個包含GPS 測站全球分布的合并單日松弛解。最后,以全球單日松弛解做為準觀測值,利用GLOBK進行卡爾曼濾波參數估計,并通過IGS站作為參考站進行赫爾默特七參數相似變換,得到ITRF08參考框架下的測站坐標和速率。

BERNESE5.2軟件解算方案:采用GPSEST程序計算基線浮動解,在此基礎上,基于最小二乘估計準則,解算獨立基線。采用QIF(Quias_Ionosphere Free)方式,忽略基線之間的相關性,逐條確定基線中的相位模糊度。然后,顧及基線之間的相關性,根據相位雙差觀測方程,利用GPSEST程序計算獲得GPS測站的坐標自由解。最后,選擇合適的框架基準,選擇ITRF框架或局部框架,采用重心基準平差的方法,利用ADDNEQ2程序將觀測網依附至特定框架系中,獲得觀測網中GPS測站的坐標約束解。在對ITRF或IGS先驗坐標作為約束時,設定閾值:N向0.01 m,E向0.01 m,U向0.03 m。根據約束閾值,從框架站列表中自動剔除重心基準下殘差超限的測站。

在GPS高精度地殼形變數據處理中,為了得到準確可靠的結果,必須對結果進行驗證。因此,采用了GAMIT10.4和BERNESE5.2兩種軟件不同算法對同一數據進行了解算,通過對兩者結果時序趨勢的一致性和離散度的比較,來判斷解算結果的質量與可靠性。

2 解算結果對比分析

2.1 基線解算結果

基線處理是GNSS地殼形變高精度數據處理中最關鍵的環節之一,其結果質量好壞直接決定著最終測站坐標、速度和地殼形變信息提取的成果。在基線解算過程中,由多臺GNSS接收機在野外通過同步觀測所采集到的觀測數據,被用來確定接收機間的基線向量及其方差-協方差陣。基線解算結果除了被用于檢驗外業觀測質量和后續的網平差外,更重要的是,因基線向量提供了點與點之間的相對位置關系,并且與解算時所采用的衛星星歷同屬于一個參考系,不受參考框架和共模誤差的影響,可確定GNSS網高精度的幾何形狀和基線變化量,對捕捉地殼微動態形變信息有著巨大的優勢。利用GAMIT10.4和BERNESE5.2兩種軟件分別解算了川滇地區2010—2014年中國地殼運動觀測網絡51個GPS連續參考站基線單日解序列,其中的兩條基線SCLH-SCXJ、SCDF-SCYX的解算結果見圖1。由圖1,兩者DX、基線長度DL離散度均在5 mm之內,DY、DZ離散度在10 mm之內。兩種軟件解算結果整體吻合度較高,時序曲線趨勢變化一致??梢姡瑑煞N軟件基線解算結果基本能夠滿足高精度地殼形變監測的要求。

圖1 基線SCLH-SCXJ、SCDF-SCYX三分量與基線長度單日解時間序列

2.2 坐標解算結果

利用GAMIT10.4和BERNESE5.2解算了川滇地區陸態網62個GPS連續參考站2010—2014共5 a的單日解坐標時間序列。其中的4個測站SCLH、SCXJ、SCDF和SCYX的解算結果見圖2。由圖2可見,兩種軟件網平差后單日解坐標時間序列離散度接近,基本都在10 mm之內,中誤差均優于3 mm。兩者坐標序列及其中誤差的整體趨勢變化一致,吻合度較高。雖BERNESE解算中誤差整體略低于GAMIT,但這并不說明中誤差小者解算精度就高,因為可能兩者解算中誤差之間存在有系統偏差。且同一軟件之間的相互比較只能說明處理方案內部的一致性。而不同軟件,不同方法的解算結果之間的比較,能更客觀地評價數據處理結果的外部一致性與可靠性。

2.3 精度評定

單日解的重復率是衡量GPS觀測質量與處理結果的重要指標之一,其值越小,說明內符合精度越高,解算質量越好?;€重復率計算公式為

(1)

式中:Rc為基線的重復性統計值;n為同一基線的總的觀測時段數;Ci為一個時段所求得的基線某一分量或邊長;σci為相應于Ci分量的中誤差;Cm為基線結果的加權平均值。

按最小二乘擬合出基線重復率和基線長度之間的線性關系,得出固定誤差和比例誤差部分,作為衡量基線精度的參考標準,見下式:

圖2 測站SCLH、SCXJ、SCDF和SCYX單日解坐標時間序列及中誤差

Rc=a+bLc.

(2)

其中:Lc為基線長度;a為常數部分,是與基線長度無關的誤差部分;b為比例部分,是與基線長度成正比的比例誤差。

圖3 基線重復率

解算川滇地區各基線,利用單天解分別計算基線重復率。選取一周內的重復率結果見圖3。其中基線重復率的固定部分單位為 mm,比例部分的基線長度L單位為m。由圖3可見,單天解的相對精度達到了10-9量級。選擇網中最長的基線SCXJ-SCMB,長度為264 428.667 m。經計算,N、E、U向和長度L的重復率分別為0.77 mm、1.14 mm、3.54 mm和1.04 mm。這說明對于最長的基線來說,其絕對重復率也能達到水平向1 mm左右、垂向3 mm左右的水平。

2.4 測站速度解算

對于長時間連續GPS觀測,為獲取GPS測站精確的速度場信息,對GPS長時間坐標時間序列進行趨勢項擬合獲取測站速度[7-8]。由于閃爍噪聲和隨機漫步噪聲均屬于冪律噪聲類型,而且原始時間序列的最優噪聲模型與PL+VW噪聲模型的極大似然值相差甚微。而且,冪律噪聲的水平可以反映時間序列的整體噪聲水平。因此,分別利用白噪聲模型和PL+VW噪聲模型進行了速度估計。

結果表明,所有測站殘差坐標時間序列的有色噪聲功率譜指數基本處于(-1,0)區間內??梢?,5 a的數據主要表現為閃爍噪聲。由于隨機漫步噪聲的長延時相關性,而川滇地區連續站坐標時間序列的時間過短,僅5 a時間,很可能還不足以準確量化隨機漫步噪聲的特征。不考慮時間序列中的有色噪聲會導致擬合速度精度的嚴重高估,顧及有色噪聲的影響,GPS測站擬合速度值變化不大,但是擬合速度中誤差將擴大約2~8倍。

3 速度場分析

選取歐亞板塊剛性穩定臺站作為測站速度計算的參考框架,解算歐亞板塊歐拉參數,將解算所得的ITRF2008框架下的速度場轉換為穩定的歐亞板塊參考框架下的相對速度場[9-10]。圖4為2010—2014年間由GPS獲得的川滇地區相對歐亞板塊地殼運動水平速度場。東西向、南北向速度估算精度優于1 mm/a。

圖4 2010—2014年相對歐亞板塊川滇地區GPS水平運動速度場

由圖4可見,除滇西南地區呈現出由西向東的增大趨勢外,川滇塊體其它區域地殼活動由北往南、由西往東呈逐漸衰減趨勢。在鮮水河斷裂帶的東北地區,測站運動速度較大,平均約為20 mm/a,運動方向為東向和東南向;滇南地區測站運動速度相對較小,平均約為7 mm/a。滇東南,運動方向為南向和東南向,滇西南運動方向為南向和南西向,在滇南中部紅河斷裂帶的南部地區,測站運動方向則表現為南向。因測站SCNC、SCSN、LUZH同位于穩定的四川盆地內部, 從而呈現出三測站運動速度較小, 約5 mm/a, 且大小相同,運動方向一致。顯然處于同一穩定塊體。由以上分析可以得出,川滇地塊運動特征主要表現為持續圍繞喜馬拉雅東結點作順時針旋轉,且仍然存在著南北向的擠壓特征。從地質動力學分析: 印度板塊與歐亞板塊的碰撞、擠壓是川滇地塊巖石層水平形變的主要驅動力[1],從而引起了青藏高原的隆升。在重力勢能作用下,造成青藏高原物質東向擠出,青藏高原物質的東向擠壓造成了川滇地塊的東移[3]。再在甘青地塊向南的擠壓力及相對穩定的華南地塊的阻擋下,青藏高原東南部物質相對歐亞板塊轉向南東方向運動,川滇地塊總體向東南方向順時針旋轉運動,區域斷裂帶呈旋走滑運動特征。由上述川滇地塊內中國地殼運動觀測網絡GPS基準站水平運動速度得到的川滇地塊整體運動特征與地質動力學分析的運動結果完全相符。

斷裂活動是地殼運動的主要形式。一般來說,絕大多數地塊內部比較穩定,構造變形主要發生在地塊之間的斷裂上。地塊之間的相互作用,主要通過其間的斷裂活動來實現[11]。為此,分別以安寧河-則木河斷裂、小江斷裂和紅河斷裂為界,選用了兩側能反映震間斷裂構造特征的GPS測站,分析了斷層兩側的運動速度差異。利用實測垂直斷層剖面平行斷層的速度分量分析了斷裂帶可能的變形寬度[12]。圖5給出了2010—2014年期間三斷裂地表GPS站平行斷層的速度剖面,利用反正切函數擬合其變化趨勢,確定各斷裂帶可能的變形寬度。其中,安寧河-則木河斷裂兩測運動速度差異較大,為11.4 mm/a,閉鎖寬度約為200 km;小江斷裂兩測運動速度差異為9.5 mm/a,閉鎖寬度約為200 km。紅河斷裂兩測運動速度差異較小,約為6.4 mm/a,變形寬度較大。因研究的斷裂帶運動特征以走滑為主,因此,利用平行斷層的速度剖面可以得到斷裂帶可能的變形寬度。從而看出,小江斷裂帶、安寧河-則木河斷裂帶處于明顯的震間閉鎖狀態。而紅河斷裂的震間期變形特征與其他斷裂帶變形特征相比,閉鎖狀態不很清晰。

圖5 2010—2014年各斷裂帶平行斷層速度剖面

4 結 論

兩種軟件不同方法的解算結果之間的比較,能更客觀地評價數據處理結果的外部一致性與可靠性,是檢驗網絡工程最終結果的一種有效手段。通過兩種軟件對2010—2014年間川滇地區中國地殼運動觀測網絡GPS基準站處理結果表明,兩種軟件網平差后單日解坐標時間序列離散度接近,基本都在10 mm之內,中誤差均優于3 mm。兩者坐標序列及其中誤差的整體趨勢變化一致,吻合度較高。由獲取的川滇地區相對歐亞板塊的GPS速度場得出,川滇地區的地殼運動除滇西南地區呈現出由西向東的增大趨勢外,川滇塊體其它區域地殼活動由北往南、由西往東呈逐漸衰減趨勢。運動方向主要表現為持續圍繞喜馬拉雅東結點作順時針旋轉,且仍然存在著南北向的擠壓特征。在鮮水河斷裂帶的東北地區,測站運動速度較大,平均約為20 mm/a,運動方向為東向和東南向;滇南地區測站運動速度相對較小,平均約為7 mm/a。滇東南,運動方向為南向和東南向,滇西南運動方向為南向和南西向。在滇南中部紅河斷裂帶的南部地區,測站運動方向則表現為南向。因測站SCNC、SCSN、LUZH,同位于穩定的四川盆地內部,從而呈現出三測站運動速度較小,約5 mm/a,且大小相同,運動方向一致,顯然處于同一穩定塊體。

以安寧河-則木河斷裂、小江斷裂和紅河斷裂為界,選用了兩側能反映震間斷裂構造特征的GPS測站,獲取了2010—2014年期間三斷裂帶地表GPS站平行斷層的速度剖面,利用反正切函數擬合其變化趨勢,確定了各斷裂帶可能的變形寬度,分析了現今川滇地區震間主要斷裂帶的活動特征。其中,安寧河-則木河斷裂兩測運動速度差異較大,為11.4 mm/a,閉鎖寬度約為200 km;小江斷裂兩測運動速度差異為9.5 mm/a,閉鎖寬度約為200 km。紅河斷裂兩測運動速度差異較小,約為6.4 mm/a,變形寬度較大。兩測相對速度越大,變形寬度越大,也就是說其閉鎖程度最高,鎖定深度最深,地震危險性也較高。因研究的斷裂帶運動特征以走滑為主。從所得的斷裂帶可能的變形寬度可以看出,小江斷裂帶、安寧河-則木河斷裂帶處于明顯的震間閉鎖狀態,很好地反映了震間斷裂帶的活動性。而紅河斷裂的震間期變形特征與其他斷裂帶變形特征相比,閉鎖狀態不很清晰??赡艿脑驗椋孩偌t河斷裂帶附近的其他斷裂系統的影響導致了斷裂帶兩測變形的特征;②跨斷層測站數量太少,統計特征不明顯。因此,需要進一步結合更多的GPS測站,包括區域站、連續站,來重點刻畫主要斷裂帶的活動特征及運動規律。

感謝中國地震局提供的陸態網觀測數據及吳偉偉提供的BERNESE解算結果。

[1] 李延興,楊國華,李智,等.中國大陸活動地塊的運動與應變狀態[J].中國科學(D輯:地球科學),2003,33(增1):65-81.

[2] 喬學軍,王琪,杜瑞林.川滇地區活動地塊現今地殼形變特征[J].地球物理學報,2004(5):806-812.

[3] 徐錫偉,聞學澤,鄭榮章,等.川滇地區活動塊體最新構造變動樣式及其動力來源[J].中國科學(D輯:地球科學),2003,33(增1):151-162.

[4] 魏文薪.川滇塊體東邊界主要斷裂帶運動特性及動力學機制研究[D].北京:中國地震局地質研究所,2012.

[5] 王閻昭,王恩寧,沈正康,等.基于GPS資料約束反演川滇地區主要斷裂現今活動速率[J].中國科學(D輯:地球科學),2008,38(5):582-597.

[6] 徐克科,伍吉倉,吳偉偉.基于GNSS時空數據的瞬態無震蠕滑信息檢測[J].地球物理學報,2015,58(7):2330-2338.

[7] DAVIS J L,BENNETT R A,WERNICKE B P.Assessment of GPS velocity accuracy for the Basin and Range Geodetic Network (BARGEN)[J].Geophysical research letters,2003,30(7).

[8] WILLIAMS S D P.The effect of coloured noise on the uncertainties of rates estimated from geodetic time series[J].Journal of Geodesy,2003,76(9-10):483-494.

[9] 徐克科,伍吉倉.利用抗差估計進行剛性板塊異常臺站探測及其運動參數求取[J].大地測量與地球動力學,2014,34(2):95-99.

[10] 徐克科,伍吉倉,呂志鵬.基于假設檢驗的剛性板塊臺站選取及全球板塊運動模型的建立[J].大地測量與地球動力學,2013,33(6):121-125.

[11] ENGLAND P, MOLNAR P. Right-lateral shear and rotation as the explanation for strike-slip faulting in eastern Tibet [J].Nature,1990, 344:140-142.

[12] SAVAGE J C,SVARC J I,PRESCOTT W H.Geodetic estimates of fault slip rates in the San Francisco Bay area[J].J Geophys Res.1999,104(B3):4995-5002.

[責任編輯:劉文霞]

Caculation and analysis of present-day crustal movement velocity fieldin Sichuan-Yunnan block

XU Keke1,WU Jicang2

(1.School of Surveying and Land Information Engineering,Henan Polytechnic University,Jiaozuo 454000,China;2.College of Surveying and Geo-Information,Tongji University, Shanghai 200092,China)

Based on GAMIT and BERNESE softwares,the GPS data of 51 base stations from 2010-2014 from Crustal Movement Observation Network of China is processed and analysed.The discrete degree of both baseline and coordinate time series are better than 10 mm.The change tendency is consistent and the coincidence degree is high.The baseline repeatability is better than 1+4×10-9Lin horizontal direction and 4+6×10-9Lin vertical direction.Crustal movement velocity field of Sichuan-Yunnan region relative to Eurasia plate is obtained.The results show that except the southwest of Yunnan which velocity increases from west to east,the crustal movement characteristics of other Sichuan-Yunnan region gradually reduces from north to south, and from west to east.Still,movement direction shows clockwise rotation around the Himalayan east nodes,and strong extrusion characteristics from north and south.Considering Anning River-Zemu River fault and Xiaojiang fault and Honghe fault as the boundary respectively, the GPS station reflecting fault characteristics on both sides are chozen. The velocity difference on both sides of fault is analyzed.By the velocity profile paralleling with the fault,the arctangent function trend is fitted and the deformation width of the fault zone is obtained.

GPS;Sichuan-Yunnan block;crustal movement;velocity field;fault zone

10.19349/j.cnki.issn1006-7949.2017.03.002

2016-02-04

國家973重點基礎研究發展計劃(2013CB733304);國家自然科學基金資助項目(41404023)

徐克科(1979-),男,副教授,博士.

P228.4

A

1006-7949(2017)03-0008-06

引用著錄:徐克科,伍吉倉.川滇地塊現今地殼運動速度場解算與分析[J].測繪工程,2017,26(3):8-13,18.

猜你喜歡
特征
抓住特征巧觀察
離散型隨機變量的分布列與數字特征
具有兩個P’維非線性不可約特征標的非可解群
月震特征及與地震的對比
如何表達“特征”
被k(2≤k≤16)整除的正整數的特征
中等數學(2019年8期)2019-11-25 01:38:14
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
詈語的文化蘊含與現代特征
新聞傳播(2018年11期)2018-08-29 08:15:24
抓住特征巧觀察
基于特征篩選的模型選擇
主站蜘蛛池模板: 国产色爱av资源综合区| 国产网站一区二区三区| 亚洲免费福利视频| 国产99在线观看| 亚洲av无码久久无遮挡| 在线毛片网站| 色亚洲成人| 国产午夜一级毛片| 911亚洲精品| 免费又黄又爽又猛大片午夜| 91麻豆国产精品91久久久| 亚洲无码视频图片| 最新精品国偷自产在线| 自拍欧美亚洲| 人妖无码第一页| 久久黄色小视频| 欧美在线精品一区二区三区| 中文字幕无码av专区久久| 国产久操视频| 2020最新国产精品视频| 一区二区日韩国产精久久| 欧美亚洲一二三区| 国产精品人莉莉成在线播放| 中文字幕永久视频| 国产精品一区二区国产主播| 成年A级毛片| 午夜啪啪网| 国产一级二级在线观看| yjizz国产在线视频网| 亚洲天堂成人在线观看| 亚洲天堂区| 制服丝袜国产精品| 亚洲综合色婷婷中文字幕| 无码国产偷倩在线播放老年人| 亚洲二区视频| 亚洲精品在线91| 无码专区国产精品第一页| 操国产美女| 日本高清免费不卡视频| 国产精品成人啪精品视频| 午夜不卡视频| 国产在线一二三区| 91精品国产无线乱码在线 | 日韩专区欧美| 亚洲毛片网站| 中文字幕日韩欧美| 国产青青操| 真人高潮娇喘嗯啊在线观看| 精品无码视频在线观看| 中文字幕资源站| 99精品国产自在现线观看| 亚洲精品自在线拍| 亚洲国产精品日韩欧美一区| 国产精品免费入口视频| 亚洲中文无码av永久伊人| 国产日韩精品一区在线不卡| 国产精品jizz在线观看软件| 日韩中文字幕亚洲无线码| 成人午夜天| 欧美一区二区啪啪| 国产女人综合久久精品视| 国内精品久久九九国产精品 | 狠狠做深爱婷婷久久一区| 国产亚洲精久久久久久无码AV| 免费国产黄线在线观看| 国产精品久久久久婷婷五月| 香蕉在线视频网站| 欧类av怡春院| 日本成人一区| 91小视频在线播放| 国产一级裸网站| 日韩精品一区二区三区中文无码| 亚洲精品第一页不卡| 亚洲无码精彩视频在线观看| 91破解版在线亚洲| 亚洲国产日韩一区| 日韩欧美国产中文| 色婷婷色丁香| 久久黄色小视频| 日本免费精品| 91无码人妻精品一区二区蜜桃| 国产成人精品18|