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

海溫觀測資料的三維時空分布和偏度特征研究

2011-09-25 03:58:58王輝贊張韌安玉柱陳奕德
海洋通報 2011年2期
關鍵詞:深度特征

王輝贊,張韌,安玉柱,陳奕德

(1.解放軍理工大學氣象學院 全軍海洋水文環境數值模擬中心, 江蘇 南京 211101; 2.中國科學院大氣物理研究所 大氣科學和地球流體力學數值模擬國家重點實驗室,北京 100029)

海溫觀測資料的三維時空分布和偏度特征研究

王輝贊1,2,張韌1,2,安玉柱1,陳奕德1

(1.解放軍理工大學氣象學院 全軍海洋水文環境數值模擬中心, 江蘇 南京 211101; 2.中國科學院大氣物理研究所 大氣科學和地球流體力學數值模擬國家重點實驗室,北京 100029)

利用包含Argo和WOD05的歷史現場觀測數據集,分析研究了全球歷史海溫觀測的三維時空分布特征,并對觀測數據的偏度特征進行研究。研究結果表明:1)海溫觀測數據在全球許多海域(特別是南大洋和北冰洋)較為稀少,大多數觀測集中在北歐沿海、北美的東西部沿海和西北太平洋區域以及其他少數區域;北半球觀測數據明顯多于南半球;2)海溫具有明顯的非對稱性分布特征,其非正態分布和有偏現象是一個全球的普遍現象,且隨著海區和深度而變化;非正態分布網格數量在任意給定的深度層次上占總網格數的比率超過80%;3)海溫觀測在各水平層次總體上具有正偏的特性,但在黑潮、灣流、南極繞極流等區域呈現出負偏的特征;具有正偏分布特征的區域網格數量占全球網格總數的比例隨著深度的增加而增加,至600 m時達到75%左右。

歷史溫度觀測;非正態分布;偏度;時空分布;非對稱性

Abstract:Based on Argo and WOD05 observations,the spatial-temporal distribution characteristics and the skewness distribution of global temperature observations are studied carefully.The results are shown as follows:1) ocean temperature data in much of the world’s ocean, especially in the Southern Ocean and the Arctic Ocean, are measured sparse and most observations are mainly concentrated in the coast of northern Europe, the west and east coasts of North America, and the Northwestern Pacific Ocean, and a few other areas; the number of observations in the Northern Hemisphere is more than that in the Southern Hemisphere obviously; 2)ocean temperature has obviously asymmetric distribution characteristics, and its nonnormal distribution is a global phenomenon, which varies with region and depth; there are more than 80% of the total 3° squares number with nonnormal distribution at any given depth; 3)the skewness of the ocean temperature observation frequency histogram is positive at different depths in all, but the skewness is negative in some regions such as Kuroshio current, Gulf stream, Antarctic Circumpolar Current and so on; The number proportion of the positive skewness distributions increased with the depth increased,which even reach 75% or so at 600 meter.

Keywords:historical temperature observations; nonnormal distribution; skewness; spatial-temporal distribution; asymmetric distribution

海洋溫度是一個描述海洋物理狀態的變量,對于研究全球海洋環流和全球氣候變化具有十分重要的意義。Argo全球觀測網建立以后,海洋溫度觀測數據逐步增加。將Argo觀測資料和歷史海溫觀測資料相結合,研究海溫的時空分布特征,具有重要的科學意義。

海洋與大氣數據統計分析處理中經常假設觀測具有正態分布的概率密度函數[1,2]。然而,地球物理系統并非一定是滿足正態分布的。比如:根據線性海浪理論,波面相對靜止水面的高度為隨機量并遵從正態分布。但理論和觀測都表明,由于非線性的影響,海浪波面高度的概率分布不是正態的。越來越多的研究表明,海表溫度、海表鹽度和海面高度等要素都呈非正態分布特征[3-6]。如,許多學者注意到赤道東西太平洋SST(Sea Surface Temperature)的非對稱性特征[3,7-9];也有學者注意到灣流系統的SST距平負偏的變化特征,并認為隨機溫度平流是造成負偏的主要原因[10]。

因此,非對稱特征或非正態特征問題是地球物理系統中存在的一個內在特征。但由于深層海洋觀測資料比較少,以往對海洋狀態場的非對稱性研究工作主要集中在海表溫度、海表高度、海面風場等表層特征研究,對于次表層或更深層的非對稱性特征研究偏少。同時,前人研究工作主要基于模式資料或再分析資料進行物理變量的非對稱性研究,而模式和觀測之間在非對稱性方面存在較大差異[3]。隨著Argo溫度的觀測前所未有的增多,需要揭示和回答的科學問題是:海溫觀測數據的三維時空分布特征如何?能否利用觀測資料揭示海洋次表層以下溫度是否存在非對稱性分布特征?偏度在不同深度層和不同地理位置存在什么分布特征呢?

1 數據預處理

本文使用的海溫觀測數據包含兩部分:一部分是來自于全球Argo計劃的浮標觀測數據[11],其時間范圍為1996年1月至2008年11月;另一部分數據來自于世界大洋觀測數據(WOD05)[12],其時間范圍從1958年到2008年。WOD05數據包含11個數據子集(OSD, CTD, MBT, XBT, SUR, APB, MRB, PFL, DRB, UOR和GLD),每種子集代表不同類型的觀測儀器。WOD05數據中的PFL子集包含72%的Argo數據,但數據更新較慢,為避免重復數據,在數據分析中不包含PFL子集數據,而用從中國Argo實時資料中心下載的Argo觀測資料替代。

WOD05數據記錄的是深度,Argo剖面觀測記錄的是壓強。為統一起見,采用Saunders[13]中的公式將Argo剖面的觀測壓強轉化為相應的深度。將原始觀測剖面離散深度數據通過垂直插值到選定的43個標準層:0, 10, 20, 30, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325, 350, 375, 400, 425, 450, 475, 500, 550, 600, 700, 800, 900, 1 000, 1 100, 1 200, 1 300, 1 400, 1 500, 1 750, 2 000, 2 500, 3 000, 2 500, 4 000, 4 500, 5000, 5500 m。由于僅對各個單獨剖面進行垂直插值處理,下文將標準層化后的這些非格點數據仍稱為觀測數據。

為方便起見,將Argo和WOD05結合的標準層剖面數據根據各自觀測的空間位置放入 3°×3°的網格區域內并進行如下質量控制:首先,如果存在兩個或多個剖面觀測數據空間位置在同一個 1°×1°的區域內,時間為同一天且采用同一類型海洋探測儀器,則認為是重復的,保留其中有效觀測數較多的一個剖面,剔除其他剖面;其次,對于每一個標準層上3°×3°區域內的數據(觀測數超過14),計算其平均值和標準差,并進行三倍標準差檢驗,如果某一數據值不在“平均值±3×標準差”的區間范圍內,則剔除該數據;第三,如果某個剖面通過以上兩步質量控制后剔除的數據數多于5個,則剔除該觀測剖面。

2 海溫觀測資料的分布狀況

標準層插值處理后的海溫樣本觀測數量在不同標準層深度上隨時間變化結果表明(圖1a):海溫從海表層到2 000 m層的逐年觀測樣本數量在1958—1970年間逐年上升,在整個70年代和80年代數量保持穩定。從90年代初到2000年觀測樣本數量減少的原因可能由于相關機構未及時將數據提交至各國海洋數據中心。而從2000年開始,由于Argo計劃的實施,從10m至2 000 m左右的觀測數據迅速增長。同時,樣本觀測總數量按月份統計的分布情況表明(圖1b):觀測樣本總數量在冬季月份相比其他季節而言稍微有所減少,但總體上各個月份分布比較均勻。這表明從觀測數據揭示的總體特征與季節性變化關系不大。

在世界海洋表面大約有5 400個3°×3°屬于海洋的小區域網格(由于地形原因,海洋區域網格數量隨著深度變深而減少),所以即使每年能收集5.4萬個全球海洋觀測數據,平均每個海洋小區域網格也只有10個觀測數據。直到大約1960年后,在整個世界大洋每年海表溫度測量的數量開始超過5.4萬個海洋觀測。此外,海溫觀測數據集中的區域分布不均,主要集中在以下少數幾個區域:北歐沿海、北美的東西部沿海和西北太平洋區域以及其他少數區域(如圖2)。同時從圖中可以看到熱帶地區的TAO觀測陣對于歷史溫度觀測數量的補充較大。在全球范圍內存在許多觀測數據稀少的地區;在海洋表層占7.2%的3°×3°海洋區域網格缺乏觀測數據;16.2%的海洋網格中觀測數量不足10個(如表1)。大多數觀測稀少地區位于南半球,南大洋和北冰洋海溫觀測數據稀少。

圖1 插值到標準層后的海洋溫度樣本觀測數量在不同層深度上分別按年份(a)和月份(b)統計的分布情況,縱坐標為對數坐標Fig.1 Number distribution statistics of ocean temperature data after interpolating to the selected standard level by year (a) and month (b).Note logarithmic scale of vertical axis.

圖2 10 m層3°×3°小區域網格內海溫歷史觀測數量空間分布圖。藍色:觀測數>10 個; 綠色:觀測數>100 個;黃色:觀測數>1000個;紅色:>10000個;白色:觀測數不足10個Fig.2 Number of historical ocean temperature observations per 3°square at 10 meter.Blue represents observation number,>10 observations; green>100 observations; yellow>1000 observations; red>10000 observations;white color means the number of observations is less than 10

表1 3°×3°網格海溫觀測數量統計aTab.1 Number of 3°squares with given number of observationsa

圖3 Argo溫度觀測總量與WOD05溫度觀測總量的百分比Fig.3 Ratio of total Argo observation number to WOD05 number

隨著 Argo計劃的實施,海溫觀測資料數量得到改善的效果如何呢?由于 Argo在不同海區的改善效果不同,將全球海洋(此處不考慮北冰洋)分為7個海區:北太平洋、北大西洋、北印度洋、南太平洋、南大西洋、南印度洋和南大洋,檢查Argo觀測數量在不同海區對WOD05觀測數量的改善效果。圖3為不同海區Argo溫度觀測與WOD05溫度觀測的百分比隨深度變化圖。圖中可以看出,Argo溫度的改善數量明顯的主要為南半球位置,在海洋中層甚至超過了WOD05數據的歷史溫度觀測數量總和。改善效果呈中間層改善多,淺層和更深層改善少的態勢,其中在1 750 m深度附近改進最大。這是由于在表層,WOD05資料數量很多,所以Argo增加改善效果有限;而Argo浮標設計觀測深度目標為2000 m,但由于各種原因,部分浮標觀測深度不足2000 m,導致該層及其以下的Argo觀測數據偏少,從而在該層改善效果不明顯。但總而言之,Argo觀測對于歷史觀測剖面資料數量進行了大大地補充,尤其對次表層觀測的改善效果十分明顯。

3 非對稱性特征分析

3.1 正態分布

正態分布檢驗用于驗證隨機變量是否服從正態分布,很多理論和方法都是基于正態分布假設提出,因此對樣本進行正態分布性質檢驗非常必要。由于標準正態分布的偏度和峰度值為0,可以用偏度和峰度值對正態分布情況進行描述。對樣本進行偏度-峰度檢驗方法[14]如下:

其中,n代表樣本數據點的個數,g1和g2分別表示偏度和峰度(其計算公式見3.2節),1σ和2σ分別表示偏度和峰度的均方差,N(0,1)表示服從均值為0、均方差為1的標準正態分布。

假設要檢驗的變量是遵從正態分布的,那么對于給定的顯著水平系數α,根據公式(1)—(4)計算則認為變量不遵從正態分布。其中,為標準正態分布中的顯著性水平為對應的臨界值。

3.2 偏度和峰度

偏度和峰度是用來衡量變量的分布密度曲線形狀的數字特征參數。偏度可以用來描述變量概率密度分布的非對稱性,峰度可以用來描述分布曲線的陡度。特別地,對于遵從正態分布的變量而言,其對應的偏度和峰度值應為0[15]。

偏度是用來衡量數據在采樣平均值附近的對稱性情況。如果其計算值為正值,表明密度分布曲線的峰點在平均值的左方,反之亦然(詳見附錄)。正態分布數據的偏度為 0,任何具有對稱性特征的數據的偏度也應該接近0。

描述分布偏度特性的一個有用指標是平均值和眾數之差[16]。對于一個中等對稱分布,有如下關于平均值、中值和眾數三者的經驗關系:平均值-眾數 ≈ 3×(平均值-中值),這也表明中值相比眾數更靠近平均值。因此, Karl Pearson提出一個簡單的衡量偏斜度的近似計算公式[17]:

式中:xd為中值。一般來說,Karl Pearson近似偏度S與偏度g1值同號,故平均值與中值之差往往與偏度g1同號,這種關系將在下文中被使用。

峰度的公式是:

對于標準正態分布,其峰度值為零。此外,帶“尖”分布的函數(即分布曲線就比較陡)其峰度大于零,而“平坦”分布的函數(即分布曲線平緩)其峰度小于零。

4 海溫的正態和偏度分布特征

4.1 正態分布特征

為了解海溫數據在每個3°×3°網格區域內的觀測特征,對每個數據觀測量超過300個的3°×3°網格區域內海洋溫度觀測進行正態分布檢驗。正態分布檢驗采用3.1節的判斷方法。以50 m和350 m層為例進行正態分布檢驗結果顯示(圖4)。圖中黑色網格區域代表該網格內的溫度觀測呈正態分布;灰色則代表非正態分布。由圖可看出,觀測數據的非正態分布特性在全球海洋中很明顯,特別是在北太平洋和北美洲的東西沿岸。

圖4 3°×3°網格區域內觀測數據在全球的正態和非正態分布特征。上、下圖分別代表50 m和350 m層。黑色:正態分布;灰色:非正態分布Fig.4 Normal and nonnormal distribution characteristics of observations per 3°square at 50 meter (a) and 350 meter(b).Black, normal distribution; gray, nonnormal distribution

圖5為正態和非正態分布區域網格數量隨深度的變化及正態分布數量所占比例的變化圖。只考慮區域網格觀測數大于 300的網格(為了使得正態分布檢驗更有效),從圖中可以看出,具有非正態分布特征的區域網格數量所占的比例具有隨著深度增加而減少的趨勢。正態分布網格數量在任何給定的深度層上占總網格數(只考慮區域網格觀測數大于300的網格)的比率不超過20%。

圖5 正態和非正態分布網格數量隨深度的變化圖。點劃線:非正態分布網格數量;點線:正態分布網格數量;實線;正態分布數占總的 3°×3°網格數量(只考慮區域網格觀測數大于300的網格)的百分比隨深度的變化(其縱坐標刻度在右側顯示)。橫坐標代表深度的對數坐標Fig.5 Normal and non-normal distribution characteristics at different depths of ocean.Dash-dot line: number of non-normal distribution; dotted line: number of normal distribution; solid line: percent of normal distribution accounted for the total number of 3° squares with more than 300 observations at corresponding depth.Note logarithmic scale of horizontal axis

4.2 偏度分布特征

由3.2介紹的理論可知,嚴格正態分布樣本的偏度為0。從4.1也可以看出,溫度觀測數據在絕大多數區域是非正態分布的。為進一步考察這個現象,通過計算各層上每個3°×3°網格區域內歷史溫度觀測數據的偏度,得到圖6所示的觀測數據分別50m和500 m偏度的全球分布特征。從圖中可以看出,在50 m層觀測數據在低緯度地區呈現負偏特征,高緯度地區呈現出正偏特征,而到500 m層負偏特征的網格數量顯著減少。值得注意的是,在黑潮、灣流和南極繞極流等區域,觀測數據較多地呈現出負偏特征。在三大洋中,偏度總體的分布特征以正偏為主,該特征越深越明顯。同時,偏度分布特征還隨著緯度而變化。在海洋表層,低緯度地區溫度偏度以負偏為主,而南半球中緯度部分地區以正偏特征為主。南北半球的溫度偏度特征基本呈對稱分布。

圖6 3°×3°網格區域內觀測數據在全球50m(上)和500m(下)層的偏度分布特征。灰色:正偏;黑色:負偏Fig.6 Skewness distribution per 3°square all over the ocean at 50m (a)and 500m (b).Gray: positive skewness;black: negative skewness

圖7 全球50m海洋所有的3°×3°網格區域內觀測數據的“平均值與中值之差”的柱狀圖。為方便比較,將正值一側的柱狀圖高曲線對應反射畫到負值一側。縱坐標為對數坐標Fig.7 Histogram of the mean-median ocean temperature difference at 50 meter for all 3° squares over the globe.Note logarithmic scale on the vertical axis.For comparison, bar heights from the positive axis are reflected onto the negative axis by a dashed line

為定量化各層偏度特征,計算柱狀圖正值部分對應面積與負值部分對應面積之差,這個數值也恰好等于所有3°×3°網格區域內的平均值之和減去中位數之和,并定義該值代表該層海洋溫度偏度總體特征。圖8為海洋溫度各層偏度總體特征隨深度變化圖。從圖中可以看出,該值幾乎在所有深度上均為正值,也就是說,海洋溫度觀測頻率直方圖總體上是正偏的。值得注意的是,不能用這個表示量進行不同深度層偏度強弱的比較,因為不同標準層對應的偏度值的幅度不同。

圖8 所有3°×3°網格區域內的海洋溫度偏度總體特征隨深度變化圖(實線)。虛線代表值為0Fig.8 Mean-median ocean temperature difference for all 3°squares at different standard levels(solid line).The dashed line denotes that the averaged value of mean-median is zero

圖9 正偏和負偏分布網格數量隨深度的變化圖。點劃線:負偏網格數量;點線:正偏網格數量;實線:正偏網格數占總的3°×3°網格數量的百分比隨深度的變化(其縱坐標刻度在右側顯示)。橫坐標代表深度的對數坐標Fig.9 Skewness distribution characteristics at different by depths of ocean.Dash-dot line: number of negative skewness; dotted line: number of positive distribution; solid line: percent of positive distribution accounted for the total number of 3°squares at corresponding depth.Note logarithmic scale of horizontal axis

圖9為正偏和負偏區域網格數量隨深度的變化及正偏分布數量所占比例的變化圖。從圖中可以看出,具有正偏分布特征的區域網格數量所占的比例隨著深度的增加而增加。在深度達到600 m時,全球海洋溫度觀測正偏網格數達到全球網格總數的75%左右。

5 小 結

本文利用Argo和WOD05歷史現場觀測數據研究了全球海溫觀測數據的時空分布和非對稱性分布特征,主要結論如下:

1、海溫觀測數據在全球許多海域(特別是南大洋和北冰洋)較為稀少;大多數觀測主要集中在北歐沿海、北美的東西部沿海和西北太平洋區域以及其他少數區域;北半球觀測數據明顯多于南半球。

2、溫度觀測數據的非正態分布和有偏現象是一個全球的普遍現象。最主要的特征包括:

1)非正態分布網格數量在任何給定深度層上占總網格數(只考慮區域網格觀測數大于300的網格)的比率超過80%;

2)從偏度水平空間分布來看,溫度觀測數據在總體上是正偏的。但在黑潮、灣流、南極繞極流等區域溫度觀測數據呈現出負偏的特征;

3)從偏度垂直層變化來看,從表層至深層,具有正偏分布特征的網格數量所占的比例隨著深度的增加逐漸增加,至600 m時達到的75%左右。

附 錄

圖A1 負偏和正偏示意圖Fig.A1 Schematic diagram of negative and positive skew

負偏(negative skew,如圖A1左側圖):偏度值< 0,左側“尾巴”長一些,分布數量集中在圖的右側,稱分布是左偏(left-skewed)的。該分布會有一些相對低的值出現。如(觀測):1, 1 000, 1 001, 1 002, 1003。

正偏(positive skew,如圖A1右側圖):偏度值>0,右側“尾巴”長一些,分布數量集中在圖的左側,稱分布是right-skewed的。該分布會有一些相對高的值出現。如(觀測):1, 2, 3, 4, 100。

[1]陳上及, 馬繼瑞.海洋數據處理分析方法及其應用 [M].北京:海洋出版社,1991:92.

[2] Emery W J, Richard E T.Data Analysis Methods in Physical Oceanography [M].Amsterdam: Elsevier, 2001:196.

[3] An S, Ham Y, Kug J, et al.El Ni?o–La Ni?a Asymmetry in the Coupled Model Intercomparison Project Simulations [J].J Clim, 2005,18:2 617–2 627.

[4] Sura P, Sardeshmukh P D.A Global View of Non-Gaussian SST Variability [J].J Phys Oceanogr, 2008, 38: 639–647.

[5] Bingham F M, Howden S D, Koblinsky C J.Sea surface salinity measurements in the historical database [J].J Geophys Res, 2002,107(C12), 8019, doi:10.1029/2000JC000767.

[6] Sura P, Gille S T.Stochastic Dynamics of Sea Surface Height Variability [J].J Phys Oceanogr, 2010, 40: 1 582–1 596.

[7] Burgers G, Stephenson D B.The “normality” of El Ni?o [J].Geophys Res Lett, 1999, 26(8): 1 027–1 030.

[8] Dong B.Asymmetry between El Ni?o and La Ni?a in a Global Coupled GCM with an Eddy-Permitting Ocean Resolution [J].J Clim, 2005,18: 3 373 – 3 387.

[9] Su J, Zhang R, Li T, et al.Skewness of subsurface ocean temperature in the equatorial Pacific based on assimilated data [J].Chin J Oceanol Limnol, 2009 ,27(3):600-606.

[10]Sura P.On non-Gaussian SST variability in the Gulf Stream and other strong currents [J].Ocean Dyn, 2010, 60(1):155-170.

[11]Roemmich D, Owens W B.The Argo Project: Global ocean observations for understanding and prediction of climate variability [J].Oceanogr, 2000, 13: 45-50.

[12]Boyer T P, Antonov J I, Garcia H E, et al.World Ocean Database 2005 [CD].Levitus S, Ed, NOAA Atlas NESDIS 60, U.S.Government Printing Office, Wash D C, 2006, 190.

[13]Saunders P M, Fofonoff N P.Conversion of pressure to depth in the ocean[J].Deep Sea Res, 1976, 23:109-111.

[14]黃嘉佑.氣象統計分析與預報方法 [M].北京: 氣象出版社,1990.

[15]White G H.Skewness, kurtosis and extreme values of Northern Hemisphere geopotential heights [J].Mon Wea Rev, 1980, 108: 1 446–1 455.

[16]Pearson K.Contributions to the mathematical theory of evolution, II: Skew variation in homogeneous material [J].Philosophical Transactions of the Royal Society of London, 1895, 186: 343-414.

[17]Stuart A, Ord J K.Kendall’s advanced theory of statistics: Distribution Theory [M].London: Oxford Univeristy Press, 1994.

Three-dimensional spatial-temporal distribution and skewness distribution characteristics based on global ocean temperature observations

WANG Hui-zan1,2, ZHANG Ren1,2, AN Yu-zhu1, CHEN Yi-de1

(1.PLA Research Center of Ocean Envrionment and Numerical Modeling, Institute of Meteorology, PLA University of Science and Technology,
Nanjing 211101, China;
2.LASG, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China)

P731.11

A

1001-6932(2011)02-0127-08

2010-11-05;

2010-11-19

中加國際科技合作項目(2008DFA22230);國家重點基礎研究發展計劃( 2007CB816005);國家自然科學基金重點項目(40730843)。

王輝贊(1983- ),男,湖南瀏陽人,博士研究生,主要從事物理海洋學研究。電子郵箱:wanghuizan@126.com。

張韌(1963- ),男,四川峨眉人,博士,教授(博導),主要從事海氣相互作用研究。電子郵箱:zren63@126.com。

猜你喜歡
深度特征
抓住特征巧觀察
深度理解一元一次方程
新型冠狀病毒及其流行病學特征認識
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
深度觀察
深度觀察
深度觀察
深度觀察
抓住特征巧觀察
主站蜘蛛池模板: 亚洲综合香蕉| 无码内射在线| 91美女视频在线| 国产美女自慰在线观看| 欧美福利在线| 成人一区专区在线观看| 人妻一区二区三区无码精品一区 | 免费人成在线观看成人片 | 一级黄色片网| 激情乱人伦| 日本精品一在线观看视频| 国产精品性| 99人妻碰碰碰久久久久禁片| 久一在线视频| 欧美狠狠干| 久久一级电影| 国产原创自拍不卡第一页| 2048国产精品原创综合在线| 亚洲视频一区| 综合色在线| 国产三级毛片| 免费一级成人毛片| 国产永久在线视频| 亚洲高清无在码在线无弹窗| 777午夜精品电影免费看| 国产精品手机视频| 国产精品美女免费视频大全| 国产玖玖玖精品视频| 成年看免费观看视频拍拍| 欧美亚洲激情| 黄色网址手机国内免费在线观看 | 特级欧美视频aaaaaa| 国产91精选在线观看| 国产尹人香蕉综合在线电影 | 国产香蕉在线视频| 国产精品久久久久鬼色| 亚洲国产欧美目韩成人综合| 日本国产精品| 嫩草影院在线观看精品视频| 婷五月综合| 国产在线一区视频| 欧美一级高清片欧美国产欧美| 久久久久中文字幕精品视频| 在线视频亚洲欧美| 欧美一区二区自偷自拍视频| 毛片久久网站小视频| 久久久久无码精品| 精品一区二区无码av| 中文字幕在线一区二区在线| 国产乱子伦手机在线| 在线高清亚洲精品二区| 亚洲综合精品香蕉久久网| 天堂网亚洲系列亚洲系列| 欧美日韩在线第一页| 国产女人18毛片水真多1| 在线观看欧美精品二区| 国产亚洲一区二区三区在线| 中文字幕调教一区二区视频| 美女无遮挡拍拍拍免费视频| 久久国产精品娇妻素人| 91国语视频| 2021国产精品自拍| 精品一区二区三区四区五区| 国产va在线观看免费| 亚洲国产高清精品线久久| 国产成人禁片在线观看| 久综合日韩| 1769国产精品视频免费观看| 中文字幕在线欧美| 五月激情婷婷综合| 亚洲日本www| 毛片网站在线看| 色妞永久免费视频| 666精品国产精品亚洲| 国产视频 第一页| 九九热精品免费视频| 国产亚洲男人的天堂在线观看 | 四虎免费视频网站| 亚洲色图综合在线| 中文字幕无码电影| 亚洲欧美综合精品久久成人网| 日韩毛片免费|