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

對流層頂高對拉薩地區溫室氣體柱濃度反演的影響*

2020-07-14 09:46:30劉丹丹黃印博孫宇松2盧興吉曹振松
物理學報 2020年13期

劉丹丹 黃印博 孫宇松2) 盧興吉 曹振松?

1) (中國科學院安徽光學精密機械研究所,中國科學院大氣光學重點實驗室,合肥 230031)

2) (中國科學技術大學研究生院,科學島分院,合肥 230026)

3) (皖西學院,電氣與光電工程學院,六安 237012)

對流層頂作為對流層和平流層之間的過渡層,對大氣痕量氣體濃度反演有非常重要的影響.理論分析對流層頂高對大氣分子含量垂直分布的影響,并結合2018年8月(6—16日)拉薩觀測數據,定量分析對流層頂高的變化對氣體柱-平均摩爾分數Xgas反演的影響.結果表明: 對流層頂高的變化改變氣體分子的垂直分布,XCO2,XCH4與對流層頂高度變化呈正相關,相關系數分別為0.998和0.780;XCO與對流層頂高度變化呈負相關,相關系數為0.994;而XH2O與對流層頂高度變化的相關性非常小.對流層頂高變化3 km,XCO2,XCH4 及 XCO 誤差范圍在 8.640%,0.035% 及 0.049% 以內.觀測期間,XH2O,XCO2,XCH4 及 XCO 日平均值分別在 3432—4287,406.1~408.2,1.673—1.720 及 0.082—0.095 ppmv 之間變化.觀測數據顯示,該地區CO2,CH4氣體柱-平均摩爾分數日變化趨勢相似,XCO2與XCH4的相關系數大部分高于0.5.觀測結果可為我國研究高原溫帶地區溫室氣體濃度的時空分布及其變化規律提供參考和第一手的直接觀測數據.

1 引 言

對流層頂是對流層與平流層之間的過渡層,熱力學上定義為垂直溫度廓線上溫度最低點所對應的層高度[1].頂高的變化與眾多因素有關,例如緯度、季節以及天氣活動等,通常對流層頂高度隨緯度增高而降低.一般來說,熱帶地區對流層頂高度約 15—20 km,中緯度地區 8—14 km;夏季對流層頂的高度高于冬季,例如青藏高原夏季對流層頂平均在 17 km[2?5].1949 年,Goody 等[6?9]在研究中發現,對流層頂的形成與 CO2,H2O,O3濃度有重要的關系,臭氧總量也與對流層頂高度有關.反之,對流層頂高度的變化也改變氣體的垂直分布,從而影響其在大氣中的總含量.

利用地基傅里葉變換光譜儀遙測溫室氣體含量時,需要利用反演算法獲得溫室氣體濃度結果.國際上地基遙測的反演算法多數是基于最優估算法或最小二乘法,例如總碳柱觀測網(total carbon column observing network,TCCON)利用標準的反演算法GFIT[10]對CO2柱濃度進行反演;大氣成分觀測網 (network for the detection of atmospheric composition change,NDACC)利用中紅外波段的反演算法SFIT獲得O3,HF,CO及CH4等氣體濃度[11];德國卡爾斯魯厄理工學院開發的PROFFIT反演算法[11]反演了柏林等地區CO2,CH4等氣體的柱濃度,其反演精度與GFIT反演算法獲得的結果相當.三種反演算法的輻射傳輸模型、后向校正存在一定的差異,例如SFIT反演算法采用FSCATM輻射傳輸模型計算模擬光譜,而PROFFIT采用KOPRA輻射傳輸模型.上述反演算法從流程上可以分為兩部分: 一是基于大氣前向模型計算大氣透過率光譜,二是通過優化反演參數使得計算值與測量值達到最小或最優,最終反演獲得待測大氣分子的濃度.在反演過程中,儀器線型函數、先驗廓線、地表壓力、太陽天頂角等因素均影響氣體濃度反演的精度[12?17].例如Wunch等[17]研究了先驗廓線對CO2等氣體濃度反演的影響,并利用機載廓線準確獲得 CO2,CO,CH4等氣體分子的校準因子;Frey等[18]和Hase 等[19]研究了儀器線型函數(instrument line shape,ILS)變化對 CO2濃度反演的影響,當 ILS振幅調制效率 (amplitude modulation efficiency,ME)改變 4%,XCO2變化 0.035%.先驗廓線作為前向模型中的重要參數,決定著反演結果的精度,而對流層頂高度的變化直接改變分子濃度的垂直分布特征,并最終影響反演結果,但目前有關對流層頂高度變化對溫室氣體濃度反演的影響尚未見報道.

拉薩市位于青藏高原中部,平均海拔3650 m,是海拔最高的城市之一,屬于高原溫帶半干旱季風氣候,全年多晴朗天氣,太陽輻射強,空氣稀薄,氣溫偏低.本文利用傅里葉變換光譜儀EM27/SUN觀測了拉薩地區溫室氣體的柱-平均摩爾分數,針對該地區對流層頂高度變化對溫室氣體濃度反演影響展開了理論及定量分析,并獲得了拉薩地區H2O,CO2,CH4及 CO 四種溫室氣體的柱濃度信息,測量結果可為高原溫帶半干旱氣候溫室氣體源與匯的研究提供數據支撐和重要參考.

2 基本原理

溫室氣體濃度反演包括垂直柱濃度(vertical column density,VCD)和柱-平均摩爾分數(Dryair Mole Fraction,DMF),反演方法是基于德國卡爾斯魯厄理工學院 (Karlsruhe Institute of Technology,KIT)氣象和氣候研究所開發的PROFFIT (PROFILE FIT) 反演算法,將最優估算法和非線性逐次迭代結合,利用Tikhonov-Phillips約束條件,在對數尺度上對溫室氣體柱濃度進行反演,反演算法主要包含前向模型和后向反演兩個部分.前向模型可以表述為

其中x是n維參數,包括未知的獨立參量,例如垂直柱濃度,u代表固定參數,例如壓力、溫度等;y是m維測量值,F是非線性模型,可線性化為下列形式:

其中K是m×n維的雅可比矩陣,計算公式K=DF/Dx,x0是線性化參考點,且y0=F(x0).利用最小二乘法,使得測量信號(ymeas)與模擬信號y的差值Dy最小:

S是反演狀態參量的協方差矩陣,采用高斯-牛頓算法解決非線性并執行多次迭代,第i+1迭代:

在大氣遙感時,由于未知量多于實測量,此時需要約束條件求解方程.在反演過程中,常采取確定的廓線或者對先驗廓線進行縮放的方法.常采用Tikhonov-Phillips約束條件,規則化方程為

其中,xa是變化參量的先驗值集合,B是規則化矩陣,g是規則化參數,在約束條件下的第i+1次迭代為

由于測量儀器是低分辨率光譜儀,因此,采用縮放先驗廓線的方式進行反演,并采用氧氣作為內部標準反演氣體柱-平均摩爾分數DMFs以減小儀誤差:

其中Columngas是氣體的柱總量,ColumnO2是氧氣的柱總量,壓強和溫度廓線均來自美國國家環境預報中心 (National Centers Environmental Prediction,NCEP)的分析數據.

3 觀測與結果分析

EM27/SUN光譜儀包含干涉儀和太陽跟蹤器(圖1(a)),探測器光譜響應范圍 0.83—2.5 μm,最大光程 1.8 cm,光譜分辨率 0.5 cm–1,可實現 CO2,CO,H2O,O2及 CH4氣體的同時觀測.EM27/SUN光譜儀觀測點位于拉薩市氣象局內(經度91.135°E,緯度 29.659°N),站點海拔高度 3.625 km(見圖1(b)).光譜記錄時間與UTC標準時間相差8 h.選擇晴朗無云的天氣進行觀測采集,光譜分辨率設置為0.5 cm–1.

圖1 (a) 觀測站點位置 (拉薩市氣象局);(b) 觀測設備 (傅里葉變換光譜儀,EM27/SUN)Fig.1.(a) Observing site (Lhasa meteorological bureau);(b) FTIR spectrometer (EM27/SUN).

3.1 數據質量控制

在反演之前,為了保證數據質量,需要對上述特殊情形進行預處理,包括干涉圖的直流分量校準、傅里葉變換、相位校準以及數據重新采樣等.引入兩個數據質量控制標準,第一個標準是最大信號強度的絕對值需大于5%,即采集光譜的最大信號強度的絕對值幅度低于5%,剔除該數據;第二個標準是干涉圖的變化不能超過10%.與此同時,為了保證獲得高精度的溫室氣體濃度,還需要記錄觀測站點的海拔高度、經緯度、溫、濕、壓等參數,以及觀測站點的實時氣象參數,例如: 地表溫度、地表壓強等.此外,觀測設備的定標誤差對反演結果也有影響,故在觀測前后需要測量儀器的線型函數消除此影響,本文利用LINEFIT軟件[19]計算H2O分子譜線獲得線型函數的調制效率和相位偏差,從而降低儀器線型函數對觀測數據的影響.

3.2 對流層頂高的變化對反演的影響

3.2.1 理論分析

對流層頂高的變化改變大氣分子含量的垂直分布,即大氣分子垂直廓線.為了分析對流層頂對痕量氣體濃度反演的影響,將氣體分子濃度南北梯度的差異及對流層頂高等參數引入到先驗廓線中,首先引入下列三個參量:

其中下標ref表示參考點,obs表示觀測點,grad_lat(jgas)是第j種氣體分子的南北半球梯度濃度,lat_ref參考廓線的緯度,lat_obs觀測點的緯度,z(i)表示大氣分層的第i層高度,ztrop是對流層頂高.第j種氣體(jgas)在第i層上(ilev)的濃度與參考廓線濃度之間的關系為

公式(10)顯示,對流層頂高對氣體廓線的調制程度取決于氣體的南北半球梯度、層高及當前層上氣體的參考廓線濃度.因此,對流層頂高對每一種氣體廓線的調制程度不同.由于50%的水汽分子集中在大約1.5 km(氣壓約850 hPa)高度以下,且 90%以 上 的 水 汽 限 制 在 5.5 km(氣 壓 約500 hPa)高度以下的大氣層中,因而對流層頂高的改變對H2O廓線幾乎沒有影響.而對集中分布在對流層附近的分子氣體廓線調制較大,例如CO2、CH4分子等.

3.2.2 定量分析

將對流層頂高由 10 km 增加到 18 km,獲得不同對流層頂高度的先驗廓線,并作為后續反演的先驗值.CO2,CO,CH4,O2以及 H2O 的濃度反演波段如表1所列.目標氣體反演波段存在一種或多種干擾分子,采用縮放干擾分子先驗廓線的方式扣除其影響.

由于目前缺少我國的大氣模式,在輻射傳輸和遙感探測領域,常用6種標準大氣模式(1976年美國標準大氣,中緯度夏季、中緯度冬季、副極地夏季、副極地冬季、熱帶)進行相關計算與反演,本文主要基于TCCON網先驗廓線,利用3.2.1節的廓線理論,獲得反演過程中涉及的8種氣體分子廓線,如圖2所示.

表1 反演波段Table 1.Inversion band.

由于分子在大氣層中的分布不同,對流層頂高的變化對不同分子的垂直分布影響不同.例如H2O,HDO分子 (圖2(a),圖2(b)是 HDO的先驗廓線乘以豐度比后的結果)主要分布在10 km以下,對流層頂高的變化對其廓線改變非常小,但CO2分子廓線受對流層頂高變化較大.從圖2可以看出,對流層頂從 10 km 增加到 18 km,CO2分子10 km 以 上 的 廓 線 上 移 (圖2(c)),CH4(圖2(d))及N2O(圖2(f))分子含量分布(廓線)隨對流層頂高變化趨勢類似,均在10—30 km高度范圍內有變化,而對流層頂高變化對CO廓線的影響主要集中在 7—20 km 高度上.通常情況下,O2分子在大氣中均勻混合,其含量是一個定值(0.2095 ppmv).分析對流層頂高的變化對溫室氣體濃度的影響,選取2018年8月16日的觀測數據進行分析,當天天氣晴朗、少云,且采集數據時間長,如圖3所示.

圖2 不同對流層頂高度下的 H2O (a),HDO (b),CO2 (c),CH4 (d),CO (e),N2O (f),HF (g) 及 O2 (h) 廓線Fig.2.The profiles of H2O (a),HDO (b),CO2 (c),CH4 (d),CO (e),N2O (f),HF (g) and O2 (h) at different heights of top troposphere.

觀測結果顯示,上午 11:30到下午 17:30,H2O氣體分子含量逐漸降低,日變化與對流層頂變化無關,如圖3(a)和圖3(e)所示.XCO2,XCH4在下午14:00達到峰值,其含量隨著對流層頂的增高而增大,且兩者的變化趨勢相似(見圖3(b)和圖3(c)),絕對偏差在中午14:00最小(見圖3(f)和圖3(g));XCO隨著時間逐漸減小,XCO與對流層頂變化的關系與CO2,CH4正好相反,其含量隨著對流層頂的增高而減小,進一步分析對流層頂變化對XCO2,XCH4及XCO的影響,將平均摩爾分數與對流層頂高進行線性擬合,如圖4所示.

XCO2和XCH4與對流層頂高變化呈正相關,XCO與對流層頂高變化呈負相關,CO2及CO日平均摩爾分數與對流層頂的相關性較高,相關系數分別為 0.998、0.994,CH4的相關性稍低,相關系數為0.78.依據線性擬合公式,對流層頂高變化1 km,XCO2變化 2.880%,XCH4變化 0.012%,XCO變化0.016%,由于一天的對流層頂高變化不超過3 km,由對流層頂高引起的XCO2,XCH4及XCO 誤差 分別為 8.640%,0.035% 及 0.049%,對流層頂高變化對不同氣體分子的影響程度不同,由于H2O主要分布在5 km以下,所以對流層頂高的變化對 H2O幾乎沒有影響,但對 CO2,CH4及CO影響較大,尤其是CO2與CO氣體分子,其原因是對流層頂高的增加,致使CO2和CO的廓線上移大于CH4廓線(如圖2).

3.3 濃度反演結果分析

綜合以上分析,在地基反演溫室氣體濃度時,需要考慮對流層頂高度的變化,以下結果均考慮當天對流層頂高,以NCEP分析數據為基礎,獲得實驗期間的溫度廓線,通過溫度廓線獲得當日對流層頂高的信息.觀測時間2018年8月,選取天氣條件 較 好 的 觀 測 日 (8 月 6,7,8,10,12,13,14,16日),10次掃描平均獲得 1組光譜,切趾采用Norton-Beer函數,共 2300 組光譜數據.由于觀測日的天氣不穩定,EM27/SUN需要在晴天薄云條件下進行觀測,每天的觀測時間段也不同,因此對日觀測結果取平均得到觀測期間的日平均值時間序列,如圖5所示.

在觀測期間,XH2O日平均值范圍在3432—4287 ppmv內變化,在 2018 年 8 月 10 日 (180810)變化較大,其標準偏差為162.8;XCO2日平均值為 406.1—408.2 ppmv,XCH4日平均值為 1.673—1.720 ppmv,但XCO2、,XCH4日變化幅度較大,尤其是7日、8日兩天,XCO2標準偏差分別為1.62和 1.35;XCH4標準偏差分別為 0.0071和0.0084;XCO 日平均值為 0.082—0.095 ppmv,日變化較為平緩,標準偏差小于0.0027.在觀測期間,XH2O,XCO2,XCH4及XCO 的 平 均 值 分 別 為3919.70,406.87,1.689 和 0.091 ppmv.與同臺設備在敦煌地區觀測的數據差異較大,尤其是CO2和CH4氣體濃度,原因可能是拉薩地區高海拔,空氣含量低,導致其 CO2,CH4濃度均偏低,再者,在敦煌地區將EM27/SUN觀測的XH2O數據與探空氣球探測的進行了對比分析,發現兩者之間相差13%,可能是由測量時間差、反演誤差、探空氣球觀測誤差等原因引起.

圖3 XH2O (a),(e),XCO2 (b),(f),XCH4 (c),(g)及 XCO (d),(h)日變化隨對流層頂高度的變化Fig.3.The diurnal variation of XH2O (a),(e),XCO2 (b),(f),XCH4 (c),(g) and XCO (d),(h) with tropopause height.

圖4 (a) CO2,(b) CH4 及 (c) CO 平均摩爾分數與對流層頂的線性擬合Fig.4.Linear fit of the average mole fraction of CO2 (a),CH4 (b) and CO (c) to the top of the troposphere.

圖5 XH2O (a),XCO2 (b),XCH4 (c),XCO (d) 的日平均序列Fig.5.Time series of XH2O (a),XCO2 (b),XCH4 (c)and XCO (d).

3.4 XCO2與XCH4相關性分析

拉薩地區XCO2,XCH4兩種氣體分子的日平均變化趨勢相似(圖5),為了分析兩者是否具有相同的源,比較了兩天中XCO2與XCH4隨時間的變化趨勢,并計算了觀測期間兩者的相關性,如圖6和圖7所示.

2018年8 月7 日在11:00~14:00(8日13:00~15:00)時段內,天空中有云,導致信號強度幅值變化大于10%,在預處理過程中將其剔除,因此圖6在該時間段出現了數據缺失.從圖6中可以看出XCO2,XCH4變化趨勢相似,隨著時間先增大后減小.將觀測期間XCO2與XCH4結果進行了相關性分析 (見圖7),其結果發現 8月 6日、7日、8日、13日、16日相關系數高于0.5,尤其是7日、8日和13日相關系數達到0.86左右,高相關系數說明 CO2與 CH4分子具有相同的源.圖7(c)、圖7(d)、圖7(f)和圖7(h)圖呈現斜“V”字型,可能在觀測期間存在源的輸送,例如中午時段的汽車尾氣,周邊地區外來源的輸送或者大氣的垂直輸送等,且拉薩地區8月份常出現反氣旋天氣現象,可能對氣體的輸送也存在一定的影響,后續將進一步研究反氣旋天氣與氣體的輸送關系.本文初步結合觀測期間拉薩地區空氣流輸送軌跡分析CO2與CH4的源,利用NOAA開發的軌跡模型HYSPLIT分析后向72 h氣體分子輸送軌跡,如圖8所示.

從氣體分子輸送軌跡來看,3類軌跡的運輸速度相當,8月6日大氣氣流的輸送軌跡主要來自從東南方向,3類軌跡的起源均來自林芝地區,途經山南地區向西北方向輸送到達拉薩;而山南市河流眾多,水資源豐富,冰川蓄積量約 10億立方米,山南市有山羊、綿羊等動物,且盛產青稞、小麥、玉米等作物,主要林木有楊、柳、落葉松等,這些均是CH4及CO2的源.8月7日和8月14日均存在垂直方向上的氣流對流,但到達拉薩觀測點已經沉降到地面,成為了局地源,8月7—16日氣體的輸送源主要來自當地的局地源(除了10—12日期間高層大氣運輸來自北邊那曲地區).

圖6 2018 年 8 月 7 日、8 日 XCO2 (a),(b),XCH4 (c),(d) 時間序列Fig.6.Time series of XCO2 (a),(b)and XCH4 (c),(d)at August 7,8,2018..

圖7 XCO2 與 XCH4 相關性 (a) 2018-08-6;(b) 2018-08-7;(c) 2018-08-8;(d) 2018-08-10;(e) 2018-08-12;(f) 2018-08-13;(g) 2018-08-14;(h) 2018-08-16Fig.7.The correlation between XCO2 and XCH4: (a) August 6,2018;(b) August 7,2018;(c) August 8,2018;(d) August 10,2018.

3.5 地基觀測結果與WACCM數據對比分析

美國大氣研究中心發展的全球大氣模型WACCM(whole atmosphere community climate model)可預測從地面到熱層頂大氣時空演化特性.在研究痕量氣體濃度方面,通常用于分析其季節變化趨勢[20],由于地基EM27/SUN的觀測時間短,在觀測期間沒有直觀的觀測數據(探空氣球、衛星數據)進行對比,只能將其與WACCM模型數據進行對比.

本文數據以拉薩為中心(經度91.135°E,緯度29.659°N),取經緯度 ± 0.5°內的數據作為拉薩地區樣本點.圖9是 WACCM模式數據與地基EM27觀測結果的對比圖,從圖9可以看出,WACCM模式獲得的拉薩地區XCO2,XCH4平均值分別為 407.76 ppmv 和 1.883 ppmv;相比較WACCM模式數據,拉薩地基觀測的XCO2值平均低 1.68 ppmv,XCH4平均低 0.193 ppmv,兩者之間的相對偏差為0.4%(XCO2)和10.24%(XCH4).推測兩者差異較大的原因主要有兩點: 一是地基觀測是通過反演算法獲得柱濃度,存在反演誤差;二是WACCM模式基于GEOS-5模型獲得的模擬數據,也會導致兩者結果具有差異.

圖8 2018 年 8 月 6 日—16 日 72 h 后向軌跡圖 (a) 2018 年 8 月 4—6 日;(b) 2018 年 8 月 7—9 日;(c) 2018 年 8 月 10—12 日;(d) 2018年 8月 14—16日Fig.8.72-hour back trajectories of Lhasa during August 6—16,2018: (a) August 4-6,2018;(b) August 7-9,2018;(c) August 10-12,2018;(d) August 14-16,2018.

圖9 地基觀測 XCO2 (a),XCH4 (b) 日平均值與 WACCM 數據對比Fig.9.Comparison of XCO2 (a) and XCH4 (b) based on ground-based observations and WACCM data.

4 結 論

基于地基便攜式傅里葉變換光譜儀EM27/SUN,獲得了拉薩地區2018年8月份(觀測時間15天)的太陽吸收光譜,反演得到 2018年8月6 日到 8 月 17 日大氣中的 H2O,CO2,CH4及 CO濃度信息.理論分析了對流層頂高對氣體分子垂直分布的影響,并定量分析了對流層頂高度變化對氣體柱平均摩爾分數的影響,分析了拉薩地區XH2O,XCO2,XCH4及XCO 的變化趨勢,并將其結果與WACCM模式數據進行了對比分析,以期為研究我國高原溫帶地區溫室氣體濃度的時空分布及其變化規律提供數據支持及理論依據,主要結果如下:

1)對流層頂高度的變化對不同氣體濃度反演的影響程度不同.H2O濃度與對流層頂高度變化沒有明顯關系;對流層頂高度變化對CO2,CH4及CO 濃度影響較大,XCO2,XCH4與對流層頂高度變化呈正相關,XCO與對流層頂高度變化呈負相關,XCO2,XCH4及XCO 與對流層頂高度變化的相關性系數分別為0.998,0.78及0.994.當對流層頂高度變化 3 km,DXCO2、DXCH4及 DXCO 分別為8.640%、0.035%及0.049%。

2)在觀測期間,拉薩地區XH2O,XCO2及XCH4日變化幅度較大,且XCO2,XCH4變化趨勢相似,多數天相關系數高于 0.5,說明 CO2,CH4具有相似的源,觀測期間空氣流主要來自山南地區的輸送;XH2O,XCO2,XCH4及XCO 日平均值分別在 3432—4287 ppmv,406.1—408.2 ppmv,1.673—1.720 ppmv以及 0.082—0.095 ppmv之間變化;XH2O,XCO2,XCH4及XCO 的 平 均 值 分 別 為3919.70,406.87,1.689 及 0.091 ppmv.

3)拉薩地區WACCM模式獲得的XCO2,XCH4與地基觀測值存在偏差,相對偏差分別為0.4%(XCO2) 和 10.24%(XCH4).相比較 WACCM 模擬值,地基觀測的XCO2,XCH4結果均偏小,地基XCO2平均低 1.68 ppmv,XCH4平 均 低 0.193 ppmv.WACCM模式數據可能需要考慮海拔高度的影響.

主站蜘蛛池模板: 亚洲中文字幕23页在线| 欧美三级日韩三级| 国产精品免费久久久久影院无码| 成人福利在线观看| 一级一级一片免费| 在线观看精品自拍视频| 日韩第九页| 女人18毛片久久| 国产交换配偶在线视频| 看av免费毛片手机播放| 在线观看av永久| 国产精品免费露脸视频| 中文字幕亚洲另类天堂| 日本午夜三级| 婷婷丁香在线观看| 中国国产A一级毛片| 午夜视频在线观看免费网站| 久久久受www免费人成| 久久天天躁狠狠躁夜夜躁| 国产在线视频二区| 欧美精品v欧洲精品| 久久国产乱子| 国产视频大全| 欧美人与牲动交a欧美精品| 国产欧美高清| 香蕉综合在线视频91| 中文字幕资源站| 中文字幕亚洲乱码熟女1区2区| 亚洲天堂视频在线观看| 不卡网亚洲无码| 亚洲天堂网2014| 老熟妇喷水一区二区三区| 亚洲色图另类| 免费国产小视频在线观看| 在线观看国产小视频| 亚洲欧美国产五月天综合| 国产精品亚洲а∨天堂免下载| 亚洲有无码中文网| 色综合网址| av午夜福利一片免费看| 无遮挡一级毛片呦女视频| 国产中文在线亚洲精品官网| 国产欧美自拍视频| 毛片网站观看| 真人高潮娇喘嗯啊在线观看| 国产h视频免费观看| 亚洲天堂视频在线免费观看| 亚洲国产成人久久精品软件 | 久久久久国产精品免费免费不卡| 欧美日韩中文国产| www.91在线播放| 亚洲高清资源| 久久久久亚洲Av片无码观看| 国产网友愉拍精品视频| 视频国产精品丝袜第一页| 久久成人免费| 国产在线精彩视频二区| 2021国产精品自产拍在线| 午夜a视频| 欧美精品在线免费| 国产精品污视频| 国产在线97| 久久99这里精品8国产| 伊人色综合久久天天| 伊人久久大线影院首页| 亚洲码一区二区三区| 久久久久无码国产精品不卡| 97se亚洲综合不卡| 国产精品理论片| 香蕉视频国产精品人| 亚洲日韩精品综合在线一区二区| 亚洲欧美在线看片AI| 狠狠亚洲婷婷综合色香| 特级欧美视频aaaaaa| 欧美精品在线视频观看| 欧美另类一区| 美女内射视频WWW网站午夜 | 国产午夜精品一区二区三区软件| 亚洲精品第1页| 美女国内精品自产拍在线播放| 国产小视频免费观看| 色视频久久|