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

基于高頻地波雷達的多時間尺度海流研究

2022-07-20 05:35:20張景思吳倫宇葛建忠丁平興
海洋學報 2022年6期
關鍵詞:風速區域

張景思,吳倫宇,葛建忠,丁平興

(1.華東師范大學 河口海岸學國家重點實驗室,上海 200241;2.國家海洋環境預報中心,北京 100081)

1 引言

海流是海洋學研究中最基礎的觀測要素之一,對海上救援以及溢油、藻類等的輸運有著十分重要的影響[1-2]。以往關于海流資料的獲取,通常采用的是船舶、浮標、潛標等傳統方法[3];海洋衛星遙感觀測通常利用熱紅外或者合成孔徑雷達進行成像,再通過對圖像的反演來獲取海流資料[4-5]。傳統觀測的優點是能夠獲取高時間分辨率的海流數據,但往往只能反映某一很小范圍的狀態[6],很難反映海流的空間變化。海洋衛星遙感的優點是能夠獲取大范圍區域的流場資料,但缺點也較為明顯,如數據易受到云層厚度影響,時間分辨率較低等[5]。而高頻地波雷達作為一種較新的觀測技術,可以利用高頻電磁波(3~30 MHz)在海表傳播時衰減小的特點以及海表對電磁波的散射機制,實現對100~400 km 范圍內的大面積流場進行觀測[7-8]。高頻地波雷達能夠獲取觀測范圍內的高時間分辨率表層海流數據,很大程度上彌補了傳統觀測手段的不足,從而兼顧了較大的空間范圍覆蓋和較高的時間采樣率,在海洋觀測中有十分重要的應用價值。目前,許多國家都在沿岸建設了雷達觀測網絡,進而實現對海流的實時監測,如美國已在其東西部沿岸區域布設了雷達觀測網絡,歐洲已部署50 多部高頻地波雷達觀測站,并在努力構建一個綜合觀測系統[9-10]。

國內外學者已經對高頻地波雷達進行了非常多的研究和應用,包括理論研究[11]、徑向流矢量合成[12]以及驗證測量結果的合理性[13]。這些研究極大地促進了高頻地波雷達的發展,同時也使高頻地波雷達廣泛地應用于數據同化[14]、溢油擴散[15]、粒子軌跡追蹤[16]等領域。這些研究表明,高頻地波雷達的發展正日趨完善,對于海洋方面的研究具有廣泛的應用價值。

高頻地波雷達長時間大范圍的觀測在海流科學研究方面具有巨大優勢,朱大勇等[17]使用一整年的高頻地波雷達資料并結合ADCP 數據,發現臺灣海峽沿岸流的季節性波動受季風控制,同時還發現了常年存在的背景流;Ebuchi 等[18]利用高頻地波雷達研究了Soya 灣暖流的季節性變化;Gough 等[19]對Farallones灣海流進行了潮周期和月際的研究,確定了該地區3 種主要的季節性環流體系。這些研究中很少有關于年際和極端事件方面的研究,這對于具有多種特征的動力過程的海流來說缺少整體性,很難全面了解海流在不同時間尺度下的動力特征和變化過程。我國近海每年都會有多個臺風登陸,對沿海區域生命財產和社會經濟造成巨大影響[20]。高頻地波雷達的優勢之一在于能夠實現全天候觀測,受天氣影響較小,目前已有許多學者對其觀測數據的準確性進行過研究[21-22],研究結果表明,高頻地波雷達除了在常規時期能夠獲取有效數據外,在極端事件時期獲取的數據也十分可靠。

近海海流受多種大氣、海洋動力過程影響,近岸復雜的地形和岸線對其影響也不可忽視,這些因素的非線性作用導致近海海流在空間、時間的不同尺度上呈現出復雜的變化。本文通過布設在朱家尖-嵊山的高頻地波雷達觀測站多年所測流場資料,對舟山外海海流在潮周期、季節、年際以及極端事件下的動力過程進行解譯,并分析海流在不同時間尺度下的動力特征和影響機制,為舟山海域赤潮、馬尾藻、漁業生產和環境監測提供基礎海流背景。

2 研究區域

舟山海域作為東海陸架的重要組成部分,地處杭州灣東側、長江口東南海域(圖1),是長江口水體南下的必經之路,同時是臺灣暖流和浙閩沿岸流交替影響的主要區域[23],由于各種水體在此相互交匯,從而形成了我國最重要的近海漁場。舟山外海地形復雜,島嶼眾多,水深變化劇烈,近岸水深約10 m,外海最深可達100 m(圖1a)。此外,舟山海域還受到季風以及臺風的影響,在這些因素的共同作用下,此處的海流動力過程復雜。

圖1 研究區域及舟山海域水深分布(a)和兩個雷達觀測站位置(b)Fig.1 The study area and the water depth in the Zhoushan sea area (a) and two radar stations (b)

3 數據和方法

3.1 高頻地波雷達資料

朱家 尖(29.893°N,122.428°E)-嵊 山(30.702°N,122.836°E)地波雷達組(圖1b)自2011 年建成以來,已累積了大量的觀測數據,本研究將以這些觀測數據為基礎,研究舟山外海域海流的多時間尺度動力過程。

本研究采用的數據來自武漢大學研發的OSMAR系統構造的真實海表流觀測場,該系統在極坐標系下利用自然三次樣條將徑向流插值到兩個扇形區域相交的公共網格上并合成得到矢量流[12]。時間分辨率為10 min,空間分辨率為0.05°,覆蓋范圍為28.95°~31.05°N,122.45°~124.45°E,觀測區域被劃分為41×43 的網格,實際觀測網格如圖1b 所示。關于該組高頻地波雷達的可靠性分析,前人已經做了較多的研究,沈志奔[8]通過對比驗證表明其在常態天氣條件下與ADCP 測流結果具有較好的一致性,其誤差在可以接受的范圍。李程等[22]對臺風“燦鴻”期間的雷達資料進行分析,對比浮標ADCP 數據發現流速平均誤差為3.43 cm/s,流向平均誤差為9.63°,相關系數分別為0.83 和0.92,具有較好的一致性。可以認定該組地波雷達的觀測數據是可信的。

本文獲取了舟山海域多個年份的雷達觀測資料(2017-2019 年),并統計各個月份的數據豐度(圖2)。其中,2017 年數據較為完整,只有7 月、8 月、10 月份豐度較低,而2018 年與2019 年均有較多月份的數據缺失。因此在后續的研究中只選取具有代表性月份的數據來研究不同季節的海流特征。

圖2 2017-2019 年各月份數據豐度百分比Fig.2 Percentage of data available for each month of 2017,2018 and 2019

3.2 風場資料

本文所使用的風場資料為歐洲中期天氣預報中心(ECMWF)下載的氣象再分析數據,選用的是距海表面10 m 高度的風速,時間分辨率為1 h,空間分辨率 為0.25°×0.25°,時間從2017年1月1日至2019年12月31日。

3.3 質量控制

本文所獲取的為已完成矢量合成的流場資料,通常由兩個及以上基站所測徑向流合成,當兩個徑向流的觀測角度夾角接近180°和0°時,會使合成的結果存在一定誤差,而當角度在90°左右時合成結果最佳[21]。此外由于天氣以及設備維護等原因,也會導致數據存在著質量問題,因此對原始數據進行質量控制是十分有必要的。本文將通過以下方法對數據進行質量控制[8]:

(1)在空間上剔除異常值,即在空間上剔除相對于平均值離散程度較大的值

式中,vt(x,y) 為t時刻空間上任意點的流速;vˉ為t時刻空間上流速的平均值;δ為標準差;N為t時刻的空間樣本容量。通過該方法定義特定時刻任意位置流速與整場平均流速的差大于兩倍標準差的流速為空間異常值[8,21]。

(2)在時間序列上剔除異常值,即在時間序列上剔除相對于平均值離散程度較大的值

式中,vx,y(t)為在(x,y)坐標下t時刻的流速;為在(x,y) 坐標下的樣本平均值;δ 為標準差;N為該坐標下的時間序列樣本容量。通過該方法定義特定位置任意時刻流速與全時段平均流速的差大于兩倍標準差的流速為時間序列異常值。

(3)在時間序列上剔除跳躍值,即剔除相鄰時間數據之間差異存在異常變幅的值

式中,vx,y(t)為在(x,y)坐標下t時刻的流速;為在(x,y) 坐標下的樣本平均值;δ 為標準差;N為該坐標下的時間序列樣本容量,該方法可以去除時間相鄰流速差異大于兩倍標準差的流速跳躍值。

以2017 年1 月的數據為例(圖3),圖中黑色矩形框為數據豐度大于60%的點位。質量控制前、后的差異主要體現在觀測區域的邊緣位置(圖3a,圖3b)。如前文提到的,這主要是合成角度造成的誤差,經過質量控制后,這些誤差得以消除。而在觀測核心區域的數據豐度(黑色矩形)所受影響并不顯著。對比質量控制前后每日的數據豐度(圖3c),發現質量控制對整體數據的豐度影響并不大,觀測區域邊緣位置數據的剔除占主要原因。可以認為質量控制對異常值的剔除效果較為明顯,對核心區域數據豐度影響不大,本文認為質量控制能夠有效剔除數據中的異常值,提升數據的可靠性。

圖3 質量控制前后數據豐度對比Fig.3 Percentage of data available before and after quality control

3.4 數據處理方法

為了確保在后續的研究過程中數據的完整性,減少觀測邊緣的合成誤差,以及數據在時間以及空間上的完整性,本文選取質量控制后數據豐度大于60%的點位進行研究。同時,本文還對數據進行了逐小時平均從而減少觀測中的隨機誤差[8]。

T_tide 是在最小二乘法的基礎上開發的一套matlab 工具包[24],能夠對潮汐潮流進行調和分析。利用T_tide 對高頻地波雷達數據進行分析,得到4 個主要分潮(M2、S2、K1、O1)的橢圓要素(最大流速、最小流速及最大流速方向等)。

余流是指海流去除潮流等周期性信號后剩下的部分,直接標示著水體的凈輸運趨勢,對海水中懸浮物和可溶性物質的稀釋、擴散及輸運等有著十分重要的指示意義。本文采用PL64 濾波器的33 h 半振幅周期低通濾波對逐小時平均的高頻地波雷達數據進行處理[25],從而濾掉半日及全日潮潮流,得到余流。

4 結果分析

4.1 潮周期變化

根據潮流調和分析得到各分潮潮流橢圓,見圖4。4 個主要分潮中,M2分潮所占比重最大,S2分潮次之,M2與S2分潮橢圓方向分布較一致(圖4a),但M2分潮橢圓大小為S2分潮橢圓的2 倍以上(注意圖中二者比例尺不同)。舟山海域以M2分潮為主導,其最大流速(橢圓長半軸)分布從東北區域向西南區域逐漸減小,最大流速方向(橢圓長軸方向)也隨著離岸距離的減小而發生偏轉,在近岸區域表現出順岸的趨勢,呈現往復流為主的特征,而外海旋轉流特征較為顯著。K1與O1分潮橢圓方向分布在外海區域基本一致(圖4b),但在近岸出現了較大的差異,近岸K1分潮振幅明顯大于O1分潮。日分潮在近岸也出現了往復流的趨勢,這可能是受到了島嶼岸線的影響。

圖4 M2、S2 潮流橢圓(a)和 K1、O1 潮流橢圓(b)Fig.4 M2,S2 (a),K1,O1 (b) tidal current ellipses

4.1.1 潮流性質

潮流的研究主要有潮流的性質以及潮流的運動形式。潮流性質是指日分潮K1、O1與半日分潮M2之間的相對比重F,可根據公式(4)判斷潮流的類型:

潮流的運動形式可以根據潮流橢圓的短軸與長軸比值K來確定,當 |K|大于0.25 時,潮流表現出較強的旋轉性,同時根據K值正負可以確定潮流的旋轉方向。由于舟山海域由M2分潮主導,因此計算M2分潮K值來判斷其運動形式。按照《港口工程技術規范(上冊》[26]中的相關規定,根據全日分潮O1、K1和半日分潮M2、S2的橢圓的長半軸可以計算最大可能流速:

圖5a 展示了舟山海域的F值分布,整個觀測區域的F值都在0.5 以下,根據《港口工程技術規范(上冊)》[26]規定,該區域屬于正規半日潮區域。觀察K值分布(圖5b)發現觀測區域內的 |K|均大于0.25,可以確定舟山外海域潮流為旋轉流,且呈順時針旋轉。根據公式(5)計算最大可能流速分布如圖5c 所示,最大可能流速介于0.5~1.1 m/s 之間,其分布從東北部到西南部逐漸減小,并在近岸地帶得到加強,顯示了島嶼邊界的影響。

圖5 F 值(a)、K 值(b)和Vmax 值(c)分布Fig.5 Distribution of F-value (a),K-value (b),and Vmax-value (b)

4.1.2 月內潮流變化

以2017 年1 月份為例,選取區域內點A(30.0°N,123.5°E)(圖1b)繪制表層海流時間序列(圖6)。流速矢量的變化過程呈現出明顯的順時針旋轉的特性(圖6a);速度與流向的變化過程呈現出明顯的周期性,能夠較好地反映大小潮的情況,大潮時最大流速可達0.8 m/s,小潮時最小流速接近0 m/s(圖6b,圖6c)。

圖6 2017 年1 月流速矢量時間序列(a)、流速時間序列(b)和流向時間序列(c)Fig.6 Time series of tidal current vector (a),time series of tidal current velocity (b),and time series of tidal current direction (c) in January,2017

4.2 季節和年際變化

對各月份數據進行質量控制和濾波處理,并計算得到月平均余流與風場。由于數據在某些月份存在一定的缺失(圖2),因此統計時選取數據較全的月份作為代表性季節進行分析。各季度月平均余流和風場結果都表現出非常強的季節變化(圖7),以1 月份為代表的冬季(圖7a 至圖7c)主要受北風影響,在離岸區域流速基本在5 cm/s 以內,而在近岸區域(靠近朱家尖島附近)存在著較明顯的從北向南的沿岸流,流速可達10~20 cm/s,從近岸往東,流速呈減小的趨勢。以6 月份為代表的夏季(圖7g 至圖7i)余流流速基本在10 cm/s 左右,分布較均勻,且大部分區域流向接近正北方向。春季(圖7d 至圖7f)和秋季(圖7j 至圖7l)屬于轉換階段,春季由北風逐漸轉變為東南風,風速減弱,余流流向逐漸向北轉換;秋季由東南風轉變為北風,風速逐漸增大,余流流向逐漸向南轉換。

對比2017 年與2018 年,冬季余流在近岸區域的變化并不顯著,均表現出明顯的自北向南的沿岸流,但在離岸區域余流變化則非常明顯,尤其在123.5°E以東和30°N 以南的區域,2018 年較2017 年流速相對增大(圖7a,圖7b)。夏季時,余流場的特征基本一致,無明顯的變化(圖7g 至圖7i)。春季時,2018 年較2017 年流速有所增強,這與風速的變化是一致的(圖7d,圖7e);秋季時,余流場具有較大的差異,2017 年的余流場方向為西北向,而2019 年的余流場方向為西南向,該變化與風場的變化存在一致性(圖7j,圖7l)。

圖7 各月平均余流與風場Fig.7 Average residual currents and wind fields for each month

為了進一步探究余流變化與風之間的關系,選取點A(圖1b)繪制余流和風的時間序列,并以24 h 為窗口滑動計算風速、風向與余流流速、流向之間的相關系數(圖8)。由于常規的風向指風的來向,與流向相反,在接下來的分析中均將風向定義為風的去向。余流流速和相關性均與風速有較強的一致性,流速最大可達0.8 m/s,同時相關系數最大可接近1。但從圖中可以發現流速的變化相對于風速變化具有一定的滯后性,主要表現為T1-T2期間當風速由增大轉為減小時,流速還保持著增長的趨勢(圖8a)。這一現象在T5-T6期間以及風速相對較小的T3-T4和T7-T8期間均有體現,而反映在相關性上就會有負相關的現象產生(圖8b 至圖8d)。風速與余流流速之間總體表現為風速較大時相關性較強,風速較小時相關性普遍較低,當風速迅速變化時,會有負相關的現象產生。說明在風速較小情況下海流以非當地風強迫導致的余流為主,如一些陸架波動、潮汐余流、密度流等。風向與流向之間也存在著類似的關系(圖8b,圖8d)。但由于風向的變化頻率較快,而余流由于存在相對滯后性,導致流向變化較為平緩,方向相關系數具有較大的波動。

選取風速大于8 m/s 時期的余流流速、流向與風速、風向(圖8 中黃色背景),計算平均值和相關系數(表1),結果顯示在風速最大期間(T1-T2和T7-T8)風速與流速相關系數分別為0.90 與0.78,方向夾角的平均值分別為48.8°和47.8°,相關系數為0.62 和0.55,相關性較強。T5-T6期間風向與流向相關系數達到了0.68,但方向夾角只有22.0°,同時風速與流速之間的相關系數只有0.48(表1)。這可能是由于此期間風速和風向均發生了較快的變化,導致Ekman 輸運并沒有充分發展并穩定。T3-T4期間,風速較小且波動較大,余流流速與風速之間的相關系數為0.55,而流向與風向之間卻出現了負相關的現象(表1)。說明相比于流速,流向對于風的響應可能更加滯后。

表1 大風時期風與余流相關性統計Table 1 Correlation statistics between wind and residual currents during strong wind

圖8 2017 年與2019 年秋季余流與風的時間序列對比和相關系數Fig.8 Time series comparison and correlation between the residual current and the wind in the autumn of 2017 and 2019

4.3 極端事件

以2017 年9 月份的臺風“泰利”(201718,TALIM)為例,高頻地波雷達探測到了完整的臺風過境數據。該臺風于2017 年9 月9 日在西北太平洋生成,并于9 月13日14 時達到了超強臺風等級。臺風從東南向雷達觀測區域接近(圖9),在臺風臨近研究區域時段中心風速達到了50 m/s,觀測區域內風速最大可達23 m/s,而后往東北方向遠離觀測區域。從空間上來看,受臺風期間極端大風影響,舟山外海的平均余流場流速整體較大,可達0.62 m/s。選擇臺風靠近觀測區和離開觀測區的兩個時刻統計其流速與流向,以及在觀測區域內的平均風速與風向。在臺風臨近前(2017 年7 月15 日0 時)平均余流流速為0.37 m/s,流向為267.7°,平均風速為14.4 m/s,風向為239.2°(圖9a),而臺風遠離時刻(2017 年7 月16 日0 時)平均余流流速為0.62 m/s,流向為233.6°,平均風速為18.2 m/s,風向為207.4°(圖9b),這兩個時段平均風向與平均流向的夾角分別為28.5°與26.2°。

圖9 臺風靠近與離開時刻余流場與風場Fig.9 The residual current and wind fields of the moment when the typhoon approach and leave

從時間序列上來看,選取點A(圖1b)繪制時間序列對比(圖10),在臺風接近的過程中,流場的潮周期運動逐漸被壓制,流速雖仍然存在著周期性變化,但流向逐漸偏向風的右側,而在臺風離開后流場又逐漸恢復成漲落潮的形態(圖10b)。在臺風作用最強烈期間(圖10 中綠色背景部分),流向基本維持在210°~330°之間,流速最大可達1 m/s 以上。余流的流向隨著風向逐漸產生逆時針旋轉,在臺風作用最強期間,平均流向為249.35°,平均風向為208.49°,其夾角為40.86°,相關系數達到了0.91(圖10b),風速與余流流速的相關系數達到了0.92,說明此時臺風對表層余流過程起主要作用,兩者之間高度相關。

圖10 2017 年9 月風速、余流流速、實際流速對比(a)和風向、余流流向、實際流向對比(b)Fig.10 Comparison of wind speed,residual current velocity,actual current velocity (a),and comparison of wind direction,residual current direction,actual direction (b) in September 2017

5 討論

5.1 不同時間尺度海流的特征

舟山海域島嶼眾多、岸線曲折,海流不僅受到長江入海徑流、潮汐、風等影響,同時還受到地形影響[27]。從整體看,地波雷達解譯的海表流動特征與前人其他形式的觀測結果吻合[27-29]。本研究通過地波雷達解譯得到了較大范圍、較高頻率的表層流場,能更好地刻畫潮流空間變化特征。

舟山外海域海流存在著較為明顯的季節變化,其主要特征為冬季向南、夏季向北的余流。由于冬季長江口徑流沿岸南下,同時冬季北風強勁[28],形成冬季沿岸較強的沿岸流;而夏季臺灣暖流較為強勁[23,30],東南向為主的風[28]形成夏季向北的余流。值得注意的是陳倩等[28]認為表層余流流速受風力影響較大,導致浙江近海表層余流流速冬強夏弱,可能是由于陳倩等[28]選取的站位更靠近岸。本文研究結果表明,冬季余流流速分布具有空間上的差異,只有在近岸區域余流流速可能出現大于夏季的情況(圖7)。由此可見高頻地波雷達數據的應用在空間上具有一定的優勢。

5.2 風與余流的關系

余流除了與徑流、潮汐、地形等相關外,與風力之間存在著較大的關系[29]。在風向、風速相對較穩定的情況下,余流對風的響應比較迅速[27],在大風時期更為顯著,但余流與風之間存在著一定的響應滯后性[31],這主要與Ekman 輸運的形成和穩定有關。Ekman[32]指出在高緯度穩定的風使海流完全發展需要經過12~24 h,并且這一時間會隨著緯度的減小而減小。而表層海流由于受到風的直接作用,發展時間會比完全發展更短。如圖8 所示,當風變化迅速時,速度相關性往往會表現出與風速一致的變化趨勢,而方向相關性則并不十分一致。由于慣性,海流相對于風存在一定滯后性,方向相關性的這一現象說明海流方向可能比流速更難改變,長期持續的風可能對余流的形成起著重要作用。

6 結論

對朱家尖-嵊山測站所測雷達表層流資料進行質量控制、潮流調和分析、低通濾波和相關性分析,可以得到以下結論。

(1)舟山海域屬于典型的正規半日潮區域,潮流運動形式以順時針旋轉流為主,最大可能流速介于0.5~1.1 m/s 之間,流速大小分布為東北區域較大,西南區域較小,在近岸得到一定的增強。

(2)2017-2019 年間舟山海域余流場的年際變化并不顯著,但存在著明顯的季節變化,主要表現為流向冬季向南,夏季向北。冬季余流整體較小,但是在沿岸處存在著非常顯著的自北向南的沿岸流,其流速可達10~20 cm/s,隨著離岸距離增大,流速逐漸減小。夏季余流整體較大,流速在10 cm/s 左右,流向基本為正北向。春秋季屬于轉換時期,流場的流向與季風的風向有關。

(3)風與余流存在著較強的相關性,尤其是在大風期間,流速與風速的相關系數在0.48~0.90 之間,流向與風向的相關系數在0.55~0.68 之間,風速對速度、方向相關性有較大影響,風速較大時相關性普遍較大,速度相關性與風速之間具有很強的一致性,而方向相關性則存在著一定的波動。

(4)在臺風期間會形成由臺風主導的余流場,在臺風作用最強烈的兩個時刻,流向夾角分別為28.5°與26.2°,流速較常態天氣下顯著增大。余流與臺風之間具有極大的相關性,臺風期間平均風向與余流平均流向夾角為40.86°,風向與余流流向的相關系數達到了0.91,風速與余流流速的相關系數達到了0.92。

猜你喜歡
風速區域
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
分割區域
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
基于GARCH的短時風速預測方法
關于四色猜想
分區域
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
GE在中國發布2.3-116低風速智能風機
主站蜘蛛池模板: 日韩精品欧美国产在线| 日韩经典精品无码一区二区| 极品私人尤物在线精品首页 | 国产精品人人做人人爽人人添| 91视频精品| 国产99视频免费精品是看6| 国产女人18水真多毛片18精品 | 国产成人亚洲无码淙合青草| 中文字幕 91| 四虎在线高清无码| 美女无遮挡免费视频网站| 国产午夜一级淫片| 五月天天天色| 一本色道久久88综合日韩精品| 国产91在线|日本| 蜜桃臀无码内射一区二区三区| 久久精品人人做人人爽电影蜜月| 91久久国产综合精品女同我| 青草国产在线视频| 国产精品嫩草影院视频| 97超碰精品成人国产| 暴力调教一区二区三区| 伊人成人在线| 中文天堂在线视频| 日本高清成本人视频一区| 欧美一级在线播放| 久久一本日韩精品中文字幕屁孩| 综合色88| 99精品一区二区免费视频| 99福利视频导航| 久久国产精品麻豆系列| 久久国产乱子伦视频无卡顿| 福利在线不卡| 91po国产在线精品免费观看| 中美日韩在线网免费毛片视频 | 国产精品30p| 日本免费一区视频| 国产成人高清亚洲一区久久| 婷婷六月综合网| 国产微拍一区二区三区四区| 亚洲午夜国产精品无卡| 波多野结衣无码中文字幕在线观看一区二区 | 欧洲欧美人成免费全部视频| 亚洲性色永久网址| 97无码免费人妻超级碰碰碰| 国产凹凸一区在线观看视频| 青青青国产精品国产精品美女| 97精品伊人久久大香线蕉| 91精品专区国产盗摄| 日韩不卡免费视频| 久久精品视频一| www.日韩三级| 亚洲an第二区国产精品| 欧美精品亚洲精品日韩专区| 香蕉视频国产精品人| 激情综合激情| 午夜日b视频| 中文毛片无遮挡播放免费| 国产另类乱子伦精品免费女| 亚洲综合香蕉| 91久久国产综合精品女同我| 男女男免费视频网站国产| 久久综合丝袜长腿丝袜| 免费啪啪网址| 国产无码网站在线观看| 日日拍夜夜操| 国产波多野结衣中文在线播放| 婷婷六月天激情| 五月婷婷精品| 国产精品lululu在线观看| 97久久超碰极品视觉盛宴| 国产美女久久久久不卡| 欧美一级99在线观看国产| 热九九精品| 免费女人18毛片a级毛片视频| 自拍偷拍欧美| 国产丝袜第一页| 亚洲性色永久网址| 伊人五月丁香综合AⅤ| 亚洲成a人片7777| 97无码免费人妻超级碰碰碰| 免费一级成人毛片|