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

利用船載走航渦動相關資料分析北黃海海域大氣湍流特征*

2014-12-02 03:50:32段自強高會旺高增祥
關鍵詞:大氣風速特征

段自強,高會旺,高增祥

(中國海洋大學海洋環境與生態教育部重點實驗室,山東 青島 266100)

海洋大氣邊界層是海氣相互作用的主要區域及海洋與大氣之間能量、物質交換的通道。大氣邊界層內流體運動的最主要特點是其湍流性,因此,許多海洋大氣試驗都對大氣湍流運動特征進行了觀測研究,如“熱帶海洋與全球大氣-海氣耦合響應試驗(TOGACOARE)”,“中尺度海氣相互作用試驗(MASEX)”,“全球大氣研究計劃大西洋熱帶試驗(GATE)”,以及“南沙群島及鄰近海區綜合考察”和“南海季風試驗(SCSMEX)”等[1-3]。

除了在南海渚碧礁和永興島近海鐵塔等固定觀測平臺進行的觀測實驗之外[4-9],中國一些學者也利用船載渦動相關系統對海洋大氣邊界層湍流特征進行了研究[10-13]。利用船載儀器進行觀測時面臨的問題是船體運動和晃動對觀 測結果 的影響[14-15]。薛宇歡等[12]、高會旺等[13]分別利用船體姿態參數、濾波等方法校正了船體運動對觀測數據的影響,并討論了利用船載渦動相關系統的定點觀測資料分析大氣湍流運動特征的可行性。與定點觀測相比,走航觀測可以在很大程度上擴展研究的時空范圍,方便地獲取不同海域、不同時間、不同環境背景下的觀測資料,但是資料處理過程中需要進一步考慮船體平移以及更加劇烈的晃動帶來的影響,因此,至今利用走航觀測資料分析海洋大氣邊界層湍流特征的研究還不多見。

本文基于東方紅2號科學考察船的船載渦動相關系統,利用2006年12月~2007年1月在北黃海海域的走航和定點觀測資料,分析了海洋大氣邊界層的湍流特征,并通過走航和定點觀測資料計算結果的對比,初步分析了利用船載走航觀測方式開展海洋大氣邊界層湍流特征觀測的可行性。

1 試驗概況與資料處理方法

1.1 觀測試驗

本次觀測試驗于2006年12月30日~2007年1月16日在中國北黃海海域進行。北黃海是位于山東半島、遼東半島以及朝鮮半島之間的半封閉陸架邊緣海,平均水深38m。本文所用資料包括船載走航觀測資料、船載定點觀測資料。為了避開船速劇烈變化對湍流脈動觀測的影響,本文僅選取船速較穩定時段下的觀測資料,共獲得120個時長為30min的走航觀測樣本。另外,觀測試驗期間在一定點站位(38°6′59″N,122°12′38″E)進行了27h的連續觀測(2007年1月13日03:00~1月14日06:00),在19個定點站位進行了時長為30min的觀測,共取得73個定點觀測樣本。每個樣本獲得的位置如圖1所示。

渦動相關系統由三維超聲風速儀(Campbell,CSAT3)組成。感應船體運動的儀器為光纖羅經(iX-sea,OCTANSⅢ),內含3個光纖陀螺和3個加速度計,可以測量船體的六自由度運動姿態以及真北方位角。船體平移速度由差分GPS(Trimble)測得。各儀器的觀測數據通過多功能數據采集器(Campbell,CR5000)采集并儲存,該系統所有儀器的采樣頻率均設為10Hz。

圖1 北黃海冬季航次觀測站位圖Fig.1 The measurement stations in the North Yellow Sea

渦動相關系統安裝平臺為東方紅2號科考船。為了減小船體對湍流觀測的影響,渦動相關系統安裝在船頭桅桿上,距海面高度為16m。三維超聲風速儀的朝向與船體中軸線一致。光纖羅經與差分GPS水平安裝在船艙內部,置于船體中心軸線上。超聲風速儀與光纖羅經的在船體中心軸線方向上的相對距離為24m,垂直方向為7m。

1.2 數據處理方法

1.2.1 數據預處理 觀測系統在觀測過程中由于受船體水平運動及晃動,以及溫度、濕度等局地環境因素的影響,觀測得到的風速、氣溫等物理要素往往相對真實數值有所偏離,因此在湍流特征計算之前需要對觀測數據進行校正[16]。數據預處理過程包括野點剔除與船體運動校正。

對于野點剔除,目前沒有統一的標準,一般采取方差檢驗的方法,即當數據點偏離序列平均值為標準差n倍以上時則認為該點為野點,并用相鄰兩點數值的線性內插值將其取代。不同學者對方差倍數n的取值不同,Vickers和 Mahrt[17]與郭建俠[18]分別采用3.5倍和4倍標準差作為檢驗標準。本文選取4倍標準差作為檢驗標準,即偏離序列平均值4倍標準差的數據被剔除,且以其前后相鄰數據內插補齊。

渦動相關系統應用于海上觀測時,船基坐標系下觀測的三維風速是真實風速與船體運動速度之和。船體運動對風速觀測的影響主要來自3個方面[14-15]:(1)船體的前后顛簸、左右搖擺以及航向變化而引起的船基坐標系的瞬時傾斜;(2)船體運動引起的船基坐標系相對于地球坐標系統的角速度;(3)船基坐標系的平移速度。

Fujitani[14]與 Edson等[15]給出了利用船體姿態參數對風速觀測值進行校正的方法:

其中:Vtrue為基于地球坐標系的真實風速矢量;Vobs為觀測得到的基于船基坐標系的風速矢量;T為將船基坐標系旋轉至地球坐標系時的變換矩陣;Ω為船基坐標系的角速度矢量;M是風速儀與船體重心之間的相對位移;S為船體姿態感應器相對于船體重心的位移;M-S為風速儀相對于船體姿態感應器的位移;Vmot是船體姿態感應器處的船體平移速度,包含船體重心的平移速度與船體角速度在感應器處引起的相對平移速度。本次觀測試驗中,光纖羅經與超聲風速儀之間的相對位移M-S為(-24,0,-7)。

1.2.2 湍流脈動量計算 計算湍流特征參數需要利用物理要素的湍流脈動值,在得到湍流脈動值之前需要對風速進行坐標旋轉,并對所有物理要素進行趨勢去除。

本文采用二次坐標旋轉方法對觀測數據進行處理,其目的是通過坐標旋轉使坐標系x軸與平均水平風方向平行,從而使平均側風速度與平均垂直風速度為零[19]。第一次坐標旋轉以z軸為中心對坐標系進行旋轉,使側風向的平均風速為零,第二次旋轉是以y軸為中心軸,在x-z平面內旋轉,從而使垂直方向風速為零。

渦動相關系統直接觀測的原始數據中包含了物理要素的平均狀態量與脈動量,為了去除物理要素的平均狀態,得到表征湍流運動的脈動量,需要對觀測數據進行趨勢去除。本文采用線性趨勢去除的方法對觀測數據進行處理[15]。

1.2.3 數據質量控制 渦動相關法作為一種直接觀測湍流特征的方法,雖然進行計算時不采用參數化的假設,但是對渦動相關系統的應用進行了前提假設,即認為觀測條件應滿足:(1)穩態(無局地變化);(2)地形平坦(無氣流的輻合輻散);(3)下墊面性質均勻一致(均一性,即物理要素濃度無水平梯度)[20-22]。

由于湍流觀測的前提假設條件往往無法滿足,觀測數據會包含有湍流非平穩性以及偏離湍流特征的無效數據,因此需要對觀測樣本進行檢驗,沒有通過檢驗的樣本即被剔除。湍流運動檢驗的方法有湍流平穩性檢驗、總體湍流特征檢驗[23]。

1.3 觀測條件

走航觀測期間風向以偏西風為主,其中西風出現的頻率為25%,西北風出現的頻率為22%,平均風速為9.4m·s-1(見圖2(a))。定點觀測期間的主導風向為南風,其出現頻率為30%,平均風速為5.9m·s-1(見圖2(b))。

圖2 觀測期間風向風速分布圖Fig.2 Wind rose diagram

大氣穩定度是表征大氣湍流發展的一個重要參數,在近地層氣象研究中常引用基于莫寧-奧布霍夫相似理論得到的無量綱高度ζ=z/L來表示。表達式中z為觀測高度;L為莫寧-奧布霍夫長度:

由圖3可以看出,冬季海洋大氣層結以單一的不穩定條件居多,這是因為對冬季海洋下墊面而言,白天和夜晚海水的溫度均高于大氣溫度,從而使大氣層結總處于不穩定狀態(z/L<0)。走航觀測期間大氣穩定度的分布較集中,頻次峰值位于-0.2處,平均值為-0.31。定點觀測期間大氣穩定度分布相對分散且更偏向于不穩定狀態,頻次峰值位于-0.3處,平均值為-0.76。

圖3 大氣穩定度頻次分布圖Fig.3 Probability density distribution of atmospheric stability

2 結果分析

2.1 功率譜

功率譜描述的是湍流脈動方差在頻率上的分布。根據Kolmogorov[25]提出的湍能串級理論,在慣性副區功率譜值與頻率之間滿足-2/3次方關系。由于船體運動具有一定的周期性,它引起的風速脈動方差會分布在與其周期相對應的頻率區間,因此風速脈動功率譜一定程度上還可以反映船體運動對觀測數據的影響,可以作為一種檢驗船體運動影響及校正效果的方法,這里以垂向風速脈動值的功率譜為例進行分析(見圖4)。

圖4 船體運動校正前后垂向風速脈動方差的譜分布Fig.4 Normalized spectra of vertical wind velocity

圖4是走航和定點觀測資料的平均功率譜,可以看出,原始譜值在0.3~1Hz的范圍內(對應自然頻率n的范圍約為0.1~0.3Hz)有異常高值,說明船體運動的影響主要在此頻率區間。對走航觀測資料而言,船體運動對譜值的影響在經過校正以后可以很好地被去除,并且在慣性副區譜值與頻率滿足-2/3次方關系,這與定點觀測結果相似。但是走航與定點觀測資料的不同之處在于,在不受船體運動影響的頻率區間,定點觀測資料的譜值在校正前后幾乎沒有變化,而對走航觀測資料而言,校正前后高頻區的譜值也未發生變化,但是低頻區的譜值有增大的趨勢。Edson等[15]在觀測研究中也發現了這種現象,并認為這可能是因為船體平移的低頻變化未被完全去除,姿態校正時的坐標旋轉使其在垂向上產生了分量。

2.2 無量綱風速標準差

根據莫寧-奧布霍夫相似理論,在近地層各物理要素無量綱標準差(σx/X*)只是穩定度參數z/L的函數,且滿足關系式σx/X*=a(-z/L)1/3。其中:σx為物理要素時間序列的標準差;X*為相應物理要素特征量的值;a為常數。

圖5是3個風速分量(u,v,w)無量綱標準差與大氣穩定度之間的關系。從圖中可以看出,3個速度分量的無量綱標準差與穩定度參數之間都近似滿足相似理論的1/3次方關系。將走航觀測資料的垂向風速分量無量綱標準差的相似函數與其它觀測試驗結果比較,發現擬合結果與 Panofsky等[26]、Hedde和 Durand[27]等在陸地、其它海域的觀測結果相近。平均風方向風速分量無量綱標準差的相似函數與馬耀明等[4]在南沙海域得到的結果相近,但是與Said[28]、Hedde和 Du-rand[27]得到的結果相比要大,比Ide[29]得到的結果小。側向風速分量無量綱標準差的相似函數擬合結果與其它觀測結果也存在差異。

與垂直方向計算結果相比,水平方向無量綱風速標準差分布的離散性相對較大,且與其它觀測試驗結果之間的差異也較大,這種情況的出現與湍流運動的產生機制以及相似理論應用的前提假設有關。在基于相似理論分析近地層大氣湍流運動特征時,對大氣湍流結構有影響的地面粗糙度、邊界層內循環、中尺度擾動以及海洋觀測環境中波浪引起的側向湍流應力等因素沒有考慮。垂向風分量受地形等其它因子的影響很小,可以很好地符合相似理論應用的前提假設,因此觀測結果與相似理論的結果非常吻合[27]。而水平方向風分量易受以上影響因素的干擾,也就是說,水平方向風速分量的湍流特征受整個邊界層特征的影響,大氣穩定度不再是支配其相似性的決定性參數,相似性不如垂向風分量的好[27]。因此,本文利用垂向無量綱風速標準差作為湍流相似性判斷的主要依據。

將走航觀測結果與定點觀測結果相比可以發現,各風速分量方差的相似函數均非常接近,這說明移動觀測平臺上得到的風速分量方差特征與定點觀測得到的結果具有可比性。

結合以上分析,以垂向風速無量綱標準差作為判斷依據可以得出,走航觀測得到的風速方差可以很好地滿足相似理論,而且計算結果與定點觀測以及前人研究工作結果均具有很好的可比性,從而可以說明走航觀測可以用于海洋大氣邊界層中湍流特征的觀測。同時,為了分析風速對走航觀測的影響,本文將觀測資料按風速大?。ㄒ?0m·s-1為界)分為兩類,并對計算結果進行了比較(見圖6)。不同風速條件下觀測結果之間沒有明顯的差異,這符合相似理論的結果,即決定無量綱標準差的主要因子是穩定度參數,而風速的影響較小。同時,還可以說明不同風速條件下船體晃動程度的不同對觀測結果的影響也不明顯,因此,在不同風速條件下的走航觀測資料均可應用于風速無量綱標準差分析。

圖5 各風速分量(u,v,w)無量綱標準差與穩定度之間的關系Fig.5 Relationship between dimensionless standard deviation of longitudinal,lateral,vertical wind velocity and stability

2.3 湍流輸送特征

其中:τ為動量通量;ρ為空氣密度;CD為表征動量輸送特征的拖曳系數;U10為10m高度處風速??梢?,拖曳系數是最不確定的因素之一,往往受風速、大氣穩定度、海面狀態等因素影響,但人們對其中的作用過程還不是十分清楚。渦動相關法作為一種通過湍流脈動值的直接測量而進一步計算海氣湍流通量的方法,近年來在陸氣和海氣通量的觀測中得到了廣泛應用,并可用于估算拖曳系數。

圖6 不同風速條件下走航觀測結果比較Fig.6 Comparison of measurements during different wind conditions

本文走航資料計算的拖曳系數平均值為103CD=1.7,這 與曲紹 厚[35]在西太平洋的計算值103CD=1.52,馬耀明等[5]在南沙海域的計算結果 103CD=1.54相比稍大,與高志球等[6]在南沙海域計算得到的非中性層結條件下的拖曳系數103CD=2.00相比偏小。整體而言,計算結果與其它海域觀測結果可比,均明顯小于陸地上的觀測結果[30]。

本次觀測試驗過程中大氣穩定度變化范圍較小且樣本數較少,無法分析拖曳系數與穩定度之間的關系,同時由于無海面狀態的相關資料,這里僅分析了拖曳系數與風速之間的關系(見圖8)。可以看出,對走航和定點觀測資料而言,在風速小于10m·s-1時,拖曳系數均有隨風速增大而變小的趨勢,而在風速大于10m·s-1的條件下拖曳系數隨風速的變化較小。以往觀測研究中,一般認為拖曳系數隨風速有增大的趨勢,但這多是對風速大于5m·s-1的觀測條件而言。低風速條件下,Vickers等[31]在最近的觀測研究也發現,拖曳系數有隨風速增大而減小的特征。與以往研究相比,本文結果中拖曳系數隨風速的變化趨勢在較高風速條件下才發生轉變,這可能與下墊面等其它觀測條件有關。

圖7 拖曳系數與風速之間的關系Fig.7 Relationship between drag coefficient and wind speed

3 結論

本文主要利用北黃海冬季航次走航和定點觀測資料,計算并分析了海洋大氣邊界層的湍流特征,包括湍流脈動方差以及湍流輸送系數等,初步分析了走航觀測用于海洋大氣邊界層湍流特征研究的可行性,主要結論如下:

(1)觀測資料中,受船體運動影響,風速脈動方差的功率譜有異常高值,經過船體運動校正后可以去除。定點觀測資料在校正前后其它頻率區間沒有明顯變化,而走航觀測資料在校正以后出現低頻譜值增大的現象,這可能是因為船體平移的低頻變化未被完全去除。

(2)依據相似理論,以及與定點觀測資料、其它研究結果的對比,對走航觀測方式的可行性進行分析。結果發現,對無量綱標準差與穩定度之間的相似函數而言,其符合相似理論的-1/3次方定律,而且與定點資料、陸地以及其它海域下墊面條件下得到的相似函數可比。對拖曳系數而言,走航和定點觀測資料計算結果的大小存在差異,但是它們與風速之間具有一致的變化趨勢與相近的擬合函數??傮w而言,走航觀測方式可以用于海洋大氣邊界層中湍流特征的觀測。

雖然走航觀測資料可以反映海洋大氣邊界層的湍流運動特征,但是其仍存在一定的不確定性,主要表現為對原始風速資料進行船體運動校正處理以后,風速脈動方差出現低頻譜值增大的現象,這仍需進一步的研究和校正。

[1]Webster P J,Lukas R.TOGA COARE:The coupled ocean-atmosphere response experiment[J].Bulletin of the American Meteorological Society,1992,73(9):1377-1416.

[2]丁一匯,李崇銀,柳艷菊,等.南海季風試驗研究 [J].氣候與環境研究,2002,7(2):202-208.

[3]丁一匯,李崇銀,何金海,等.南海季風試驗與東亞夏季風 [J].氣象學報,2005,62(5):561-586.

[4]馬耀明,王介民,劉巍,等.南沙海域近海層大氣湍流結構及輸送特征研究 [J].大氣科學,1997(03):79-85.

[5]馬耀明,王介民,張慶榮,等.南沙海域大氣湍流通量輸送特征分析 [J].高原氣象,1997(01):48-52.

[6]高志球,馬耀明,王介民,等.南沙群島海域近海面粗糙度,中性拖曳系數及總體交換系數研究 [J].熱帶海洋,2000,19(1):38-43.

[7]閆俊岳,姚華棟,李江龍,等.1998年南海季風爆發期間近海面層大氣湍流結構和通量輸送的觀測研究 [J].氣候與環境研究,2000,5(4):447-458.

[8]閆俊岳,姚華棟,李江龍,等.2000年南海季風爆發前后西沙海域?!獨鉄崃拷粨Q特征 [J].海洋學報,2003,25(4):18-28.

[9]閆俊岳,唐志毅,姚華棟,等.2002年南海季風爆發前后西沙海區海-氣通量交換及其變化 [J].地球物理學報,2005,48(5):1000-1010.

[10]曲紹厚,林錫貴.西太平洋熱帶海域強對流天氣過程湍流通量輸送的某些特征 [J].高原氣象,1994,13(1):75-83.

[11]曲紹厚,胡非,李亞秋.1998年SCSMEX期間南海夏季風海氣交換的主要特征 [J].氣候與環境研究,2000,5(4):434-446.

[12]薛宇歡,高會旺,高增祥.船載渦動相關觀測系統及其有效性研究 [J].中國海洋大學學報:自然科學版,2009,39(Sup):177-183.

[13]高會旺,顧明,王仁磊,等.北黃海海域大氣湍流強度特征及風速標準差相似性分析 [J].中國海洋大學學報:自然科學版,2009,39(4):563-568.

[14]Fujitani T.Method of turbulent flux measurement on a ship by using a stable platform system [J].Pap Meteor Geophys,1985,36(3):157-170.

[15]Edson J B,Hinton A A,Prada K E,et al.Direct covariance flux estimates from mobile platforms at sea[J].Journal of Atmospheric and Oceanic Technology,1998,15:547-562.

[16]徐自為,劉紹民,宮麗娟,等.渦動相關儀觀測數據的處理與質量評價研究 [J].地球科學進展,2008,23(4):357-370.

[17]Vickers D,Mahrt L.Quality control and flux sampling problems for tower and aircraft data[J].Journal of Atmospheric and Oceanic Technology,1997,14(3):512-526.

[18]郭建俠.華北玉米下墊面湍流輸送特征及參數化方案比較研究[D].南京:南京信息工程大學,2006.

[19]Tanner C B,Thurtell G W.Anemoclinometer Measurements of Reynolds Stress and Heat Transport in the Atmospheric Surface Layer[R].Wisconsin:University of Wsconsinmedison,1969.

[20]Massman W J,Lee X.Eddy covariance flux corrections and uncertainties in long-term studies of carbon and energy exchanges[J].2002,113(1-4):121-144.

[21]Baldocchi D D.Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems:past,present and future[J].Global Change Biology,2003,9(4):479-492.

[22]Lee X,Massman W J,Law B E.Handbook of micrometeorology:aguide for surface flux measurement and analysis[M].[S.I.]:Kluwer Academic Publishers,2004.

[23]Foken T,Wichura B.Tools for quality assessment of surfacebased flux measurements[J].Agricultural and Forest Meteorology,1996,78(1-2):83-105.

[24]Fairall C W,Bradley E F,Hare J E,et al,Bulk parameterization of air-sea fluxes:Updates and verification for the COARE algorithm [J].Journal of Climate,2003,16:571-591.

[25]Kolmogorov A N.Dissipation of energy in locally isotropic turbulence[C].[S.I.]:Dokl Akad Nauk SSSR,1941:15-17.

[26]Panofsky H A,Tennekes H,Lenschow D H,et al.The characteristics of turbulent velocity components in the surface layer under convective conditions [J].Boundary-Layer Meteorology,1977,11(3):355-361.

[27]Hedde T,Durand P.Turbulence intensities and bulk coefficients in the surface layer above the sea[J].Boundary-Layer Meteorology,1994,71(4):415-432.

[28]Said.Experimental study of the marine boundary Layer:turbulence strucre and surface fluxes(TOSCANE-T Experiment)in Frence[D].Toulous France:University of Paul Sabatier,1988.

[29]Ide.dynamics and turbulence in the Sahelian surface atmospheric boundary layer(STARS experiement)[D].Toulous France:U-niversity of Paul Sabatier,1981.

[30]肖霞.黃土高原半干旱區荒草地湍流和湍流能量傳輸特征及能量平衡狀況 [D].蘭州:蘭州大學,2011.

[31]Vickers D,Mahrt L,Andreas E L,Estimates of the 10-m neutral sea surface drag coefficient from aircraft Eddy-Covariance meas-urements[J].Journal of Physical Oceanography,2013,43(2):301-310.

猜你喜歡
大氣風速特征
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
大氣古樸揮灑自如
大氣、水之后,土十條來了
新農業(2016年18期)2016-08-16 03:28:27
基于GARCH的短時風速預測方法
考慮風速分布與日非平穩性的風速數據預處理方法研究
主站蜘蛛池模板: 亚洲精品在线观看91| 国产凹凸视频在线观看| www.精品国产| 日韩欧美91| 黄色国产在线| 亚洲欧洲天堂色AV| 国产成人精品亚洲日本对白优播| AV天堂资源福利在线观看| 亚洲精品大秀视频| 免费一看一级毛片| 成人小视频在线观看免费| 国产福利免费视频| 欧美 国产 人人视频| 日韩在线网址| 欧美激情视频二区| 992Tv视频国产精品| 国产成人永久免费视频| 欧美区一区| 亚洲aⅴ天堂| 天天操天天噜| 91麻豆久久久| 亚洲精品在线91| 毛片免费高清免费| 毛片网站观看| 亚洲国产成人自拍| 五月婷婷丁香综合| 日本人又色又爽的视频| 欧美视频在线不卡| 国产高清不卡| 欧美三级不卡在线观看视频| 久久精品只有这里有| 97青青青国产在线播放| 日韩一区精品视频一区二区| 国产成人一区免费观看| 久久免费视频播放| 黄色不卡视频| 91小视频版在线观看www| 国产人碰人摸人爱免费视频| 国产草草影院18成年视频| 国产精品视频3p| 免费一级毛片完整版在线看| 国产日韩欧美在线播放| av天堂最新版在线| 日本在线欧美在线| 亚洲一区二区在线无码| 日本影院一区| 中国国产一级毛片| 精品国产aⅴ一区二区三区| 国产区成人精品视频| 四虎精品国产永久在线观看| 91久久偷偷做嫩草影院电| 欧美特黄一级大黄录像| 亚洲狠狠婷婷综合久久久久| 亚洲国产精品美女| 久草网视频在线| 色有码无码视频| 久精品色妇丰满人妻| 亚洲日韩精品无码专区| 免费无码在线观看| 国产在线观看第二页| 99热国产这里只有精品9九| 无码又爽又刺激的高潮视频| 亚洲香蕉在线| 五月天婷婷网亚洲综合在线| 亚洲人妖在线| 亚洲色图综合在线| 99热最新在线| 少妇人妻无码首页| 欧美a级完整在线观看| 亚洲中文字幕在线一区播放| 欧美日韩在线亚洲国产人| 四虎永久免费地址| 亚洲国产天堂久久综合226114| 国产99精品久久| 一级成人a毛片免费播放| 亚洲天堂啪啪| 国产91视频观看| 日韩欧美高清视频| 亚洲天堂啪啪| 欧美性精品| 免费高清自慰一区二区三区| 亚洲 欧美 偷自乱 图片|