亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        華北地區(qū)一次大氣污染過程的PM2.5/PM10同化和模擬試驗

        2017-11-07 07:48:27潘曉濱胡富嘉
        江西農(nóng)業(yè)學報 2017年11期
        關(guān)鍵詞:氣溶膠增量觀測

        黃 然,潘曉濱,李 毅,胡富嘉

        (1.國防科學技術(shù)大學 氣象海洋學院,江蘇 南京 211101;2.解放軍61867部隊,北京 101104;3.解放軍95809部隊,河北 滄州 061736)

        華北地區(qū)一次大氣污染過程的PM2.5/PM10同化和模擬試驗

        黃 然1,2,潘曉濱1,李 毅1,胡富嘉3

        (1.國防科學技術(shù)大學 氣象海洋學院,江蘇 南京 211101;2.解放軍61867部隊,北京 101104;3.解放軍95809部隊,河北 滄州 061736)

        利用WRF-Chem模式MOSAIC氣溶膠變量建立起的PM2.5和PM10氣溶膠同化系統(tǒng),對2015年12月華北地區(qū)的一次污染過程進行了資料同化和模擬試驗。結(jié)果表明,同化能夠有效改進模式初始場PM2.5和PM10的空間分布,其均方根誤差分別降低了43.7 μg/m3和71.0 μg/m3,相關(guān)系數(shù)均提高了0.52,重污染區(qū)的初始位置和強度與實況更接近。同化對于預(yù)報也有明顯的改進效應(yīng),24 h內(nèi)PM2.5和PM10模擬的均方根誤差平均降低了17.9 μg/m3和24.0 μg/m3,相關(guān)系數(shù)分別平均提高了0.20和0.18;由于排放源的持續(xù)強迫作用,初始場同化的改進效果在18 h后趨于減弱。

        WRF-Chem模式;三維變分;氣溶膠;資料同化

        0 引言

        大氣氣溶膠是懸浮在大氣中的固態(tài)或液態(tài)顆粒物和氣體載體共同組成的多項體系,粒徑在1×10-3~100 μm[1],通常劃分為PM10(<10 μm)、PM2.5(<2.5 μm)和PM1(<1 μm)。PM2.5和PM10在大氣中壽命相對較長,對氣候變化、環(huán)境污染、人體健康和農(nóng)業(yè)生產(chǎn)都有直接或間接的影響[2]。利用大氣化學模式對PM2.5和PM10進行模擬和預(yù)報是當前大氣化學領(lǐng)域的前沿研究方向之一,NAQPMS、CMAQ、WRF-Chem和GEOS-CHEM等區(qū)域或全球大氣化學模式已在氣溶膠污染數(shù)值預(yù)報領(lǐng)域得到了廣泛應(yīng)用[3]。

        盡管大氣化學模式在近幾年取得了長足進展,但由于排放源、物理化學過程、模式初始場等方面的不確定性,目前對于PM2.5和PM10的定量預(yù)報準確率仍有待提高。準確的初始場是模式預(yù)報成功的重要條件之一,資料同化(Data Assimilation)能夠?qū)⒂^測信息和模式背景場融合,給出大氣真實狀態(tài)的最優(yōu)估計[4],從而有效提高模式的預(yù)報準確率。資料同化在氣象業(yè)務(wù)上已得到普遍運用,但大氣化學模式的資料同化,特別是針對氣象化學耦合條件下的氣溶膠資料同化研究還較為滯后。

        在早期的氣溶膠資料同化研究中,同化系統(tǒng)僅僅考慮觀測資料(PM2.5或PM10)中氣溶膠總量的一種,如崔應(yīng)杰等[6]基于MM5和NAQM耦合模式,采用最優(yōu)插值方法(OI)對PM10總量、NO2和SO2分別進行資料同化,結(jié)果表明同化能較好地修正模式預(yù)報場;Denby等[7]基于LOTOS-EUROS大氣化學模式,分別采用集合卡曼濾波(EnKF)和統(tǒng)計差值(SI)兩種方法對PM10總量進行同化,并且對兩種方法的預(yù)報結(jié)果進行了比較,結(jié)果顯示SI方法對于PM10的同化更為有效。隨著氣溶膠觀測資料的豐富和氣溶膠方案的發(fā)展,在同化中按照化學成分和粒徑大小對氣溶膠進行劃分已成為主流趨勢,如Jiang等[8]基于WRF-Chem模式中的GOCART方案進行了PM10的同化試驗,將氣溶膠細化為14個變量,結(jié)果顯示12 h內(nèi)同化能有效提高預(yù)報質(zhì)量;Li等[9]基于WRF-Chem模式中的MOSAIC氣溶膠方案,將PM2.5劃分為5個物種,利用3D-Var方法對地面PM2.5進行了同化,結(jié)果表明同化對初始場和24 h內(nèi)的預(yù)報均有改進效果。

        華北地區(qū)是我國霧霾高發(fā)區(qū)域之一,其中由農(nóng)作物焚燒和工業(yè)排放等引起的空氣重污染事件極具代表意義,在冬季具有季節(jié)高發(fā)性。目前,國內(nèi)的氣溶膠觀測資料主要是PM2.5和PM10,而針對這兩類觀測資料的同化研究還比較少。本文利用華北地區(qū)的氣象和氣溶膠觀測資料,基于WRF-Chem模式發(fā)展建立的PM2.5或PM10氣溶膠同化系統(tǒng)[8],對2015年12月一次大氣污染過程進行了PM2.5和PM10兩類資料的同化試驗和模擬預(yù)報,旨在檢驗同化對預(yù)報的改進效果。

        1 模式和同化方案設(shè)計

        WRF-Chem是美國最新發(fā)展的區(qū)域大氣動力—化學耦合模式,是在NCAR開發(fā)的中尺度數(shù)值預(yù)報氣象模式中加入化學模塊集成而成,可同步計算氣象要素與大氣化學成分,實現(xiàn)“在線”耦合[10-15]。表1為本文研究所選用的主要物理和化學參數(shù)方案,其中氣溶膠方案采用的是由Zaveri等[16]提出的MOSAIC方案,其綜合考慮了氣溶膠的各個類型、成分(有機、無機氣溶膠及揮發(fā)性有機物等),是WRF-Chem模式中一種主流的氣溶膠方案。MOSAIC方案按照成分將氣溶膠分為8類,分別為EC(元素碳)、BC(黑碳)、NO3(硝酸鹽)、SO4(硫酸鹽)、Cl(氯化物)、Na(鈉鹽)、NH4(銨鹽)和OIN(無機物),并按照氣溶膠粒徑大小對每一類劃分為4組,依次為0.039~0.1、0.1~1.0、1.0~2.5、2.5~10 μm,前3組粒徑段之和為PM2.5,后1組粒徑段代表PM2.5~10。

        表1 主要物理/化學參數(shù)方案

        本文所采用的模式為WRF-Chem V3.5.1,模式設(shè)計為兩重嵌套區(qū)域(圖1):第一重網(wǎng)格嵌套以(109.4° E,36° N)為中心,緯向網(wǎng)格數(shù)為92,經(jīng)向網(wǎng)格數(shù)為84,網(wǎng)格距為45 km,包含了大部分中國地區(qū);第二重網(wǎng)格嵌套以(114.932° E,38.895° N)為中心,緯向網(wǎng)格數(shù)為97,經(jīng)向網(wǎng)格數(shù)為85,網(wǎng)格距為15 km,主要覆蓋華北地區(qū)。兩重嵌套網(wǎng)格在垂直方向上均分為30層,分辨率自下而上隨高度增加而逐漸降低。模式將采用NCEP提供的1°×1°、每6小時一次的FNL再分析資料作為氣象初始場和邊界條件?;瘜W初始場由前一天的模擬結(jié)果導(dǎo)入,排放源文件由清華大學制作提供(http://www.meicmodel.org/index.html)。

        圖1 WRF-Chem模擬區(qū)域

        每日逐時的氣溶膠觀測資料來自全國城市空氣質(zhì)量實時發(fā)布平臺(http://113.108.142.147;20035/emcpublish/),該平臺主要提供全國范圍內(nèi)1500余個站點的空氣質(zhì)量實時監(jiān)測數(shù)據(jù),資料包括了PM2.5、PM10、SO2、NO2、CO、O3的逐小時濃度值和24 h平均值。圖2為全國監(jiān)測站點地理分布和區(qū)域海拔高度,可從圖中看出站點分布呈現(xiàn)東密西疏分布不均的特點。

        針對不同的觀測資料、空氣質(zhì)量模式和模式氣溶膠方案,需要設(shè)計不同的氣溶膠資料同化系統(tǒng)。采用基于Li等[8]開發(fā)的面向MOSAIC氣溶膠方案發(fā)展的三維變分資料同化系統(tǒng)。同化過程中,目標泛函為:

        (1)

        公式(1)中,x是同化后的分析場,xb是背景場,通常由上一次模式預(yù)報結(jié)果導(dǎo)入;y代表觀測場;H為線性觀測算子,表示將模式網(wǎng)格的值插值到觀測位置;R代表觀測誤差協(xié)方差,B代表背景誤差協(xié)方差。這里采用的控制變量為PM2.5和PM2.5~10,即將模式變量的前3個粒徑段合并為PM2.5,模式變量的最后一個粒徑段為PM2.5~10。本文采用NMC方法,利用模式一個月的模擬結(jié)果,以每天00UTC的48 h預(yù)報場和24 h預(yù)報場作差,用差值場來統(tǒng)計背景誤差協(xié)方差。

        圖2 全國監(jiān)測站點分布及海拔高度

        為便于求解,將上式轉(zhuǎn)化為增量形式:

        (2)

        公式(2)中,δx=x-xb,代表增量場,d=y-Hxb,代表更新向量。我們求解得到同化后的分析增量,與背景場結(jié)合便得到需要的同化分析場x=xb+δx,這種形式構(gòu)造的代價函數(shù)綜合了背景場和觀測信息,得出最終所需的最優(yōu)同化分析場。

        根據(jù)上述同化方法,設(shè)計以下2種試驗方案,對華北地區(qū)一次污染過程中的PM2.5和PM10進行個例同化試驗。模擬試驗均采用相同的氣象條件,對不同化學初始場預(yù)報后進行比對,結(jié)果的討論分析主要針對D02區(qū)域進行;模擬時間為2015年12月29日00UTC至12月30日00UTC。

        方案1:控制試驗,不進行同化,以背景場作為初始場進行24 h預(yù)報,記為Control。

        方案2:同化試驗,同化區(qū)域為兩重嵌套網(wǎng)格區(qū)域,以分析場作為初始場進行24 h模擬預(yù)報,記為Da。

        2 同化和模擬預(yù)報結(jié)果分析

        2.1 初始場同化效果分析

        圖3(a)為地面PM2.5觀測值,可以看出重污染區(qū)域主要集中在河北中部-山西中南部-陜西中部一線,另山東北部也有局地的重污染區(qū)。圖3(b)為初始背景場,對于高污染區(qū)域河北-山西-陜西這一線模擬明顯偏弱,山東北部的大值區(qū)域也未體現(xiàn)。但總體而言,模式仍能較為準確地描述污染物的空間分布。圖3(c)為PM2.5同化增量場,在河北中部、山西中部和陜西中部分別有一個較強的正增量中心,增量值均超過了150 μg/m3,在山東北部和山西南部還存在相對較弱的PM2.5正增量中心,與觀測的重污染區(qū)相對應(yīng);此外,在河南北部和安徽北部有兩個較強的負增量中心,數(shù)值均在-50 μg/m3左右,與這兩個區(qū)域的低濃度相對應(yīng)。圖3(d)為同化后分析場,可以看出污染中心轉(zhuǎn)移至河北中部,濃度最大值為326 μg/m3,污染范圍增大,同時山東北部也出現(xiàn)重污染區(qū),與實況分布更為接近。

        圖4(a)為地面PM10觀測值,分布基本與PM2.5一致,在山西南部有觀測最大值693 μg/m3。圖4(b)為PM10初始背景場,可以看出唐山區(qū)域模擬值比觀測值偏高,其余區(qū)域較觀測值偏低。圖4(c)為PM10的同化分析增量,在山西中部和陜西中部有兩個正增量中心,正增量值均超過了220 μg/m3;在河北中部、山西南部和河南北部也分別出現(xiàn)了較強的正增量中心;在河北唐山為負增量場,最大值達到了-108 μg/m3。圖4(d)為PM10同化后分析場,可以看出污染中心由唐山南移至河北中部,中心濃度值達到402 μg/m3,唐山區(qū)域濃度值則降至340 μg/m3以下。

        為了定量分析同化對初始場的改進效果,采用均方根誤差(RMSE)和相關(guān)系數(shù)(CORR)對PM2.5和PM10進行檢驗分析,將分析場插值到觀測站點與觀測值進行比對,以此來驗證同化前后效果。統(tǒng)計指標分別為:

        圖5為D02區(qū)域同化前后初始場的散點圖。圖5(a)為PM2.5散點圖,從圖中可看出未經(jīng)同化時(藍色散點),散點主要分布于中心軸線兩側(cè),以右偏的情況居多,說明同化前模式背景場的PM2.5普遍偏低,特別是對于觀測值達到400 μg/m3左右的觀測點,其對應(yīng)的模式值均未超過200 μg/m3。同化后,分析場與觀測場的PM2.5(紅色散點)均密集于中心軸線兩側(cè),說明同化后分析場的PM2.5更接近觀測實況。從統(tǒng)計結(jié)果來看,RMSE由84.5 μg/m3降低至40.8 μg/m3,降低了43.7 μg/m3;CORR由0.31提高至0.83,提高了0.52。圖5(b)為PM10散點圖,散點分布情況基本與PM2.5一致,同化后的分析場比背景場更接近觀測場,RMSE由原來的125 μg/m3降低到54 μg/m3,降低了71 μg/m3;CORR由同化前的0.35提高至0.87,提高了0.52;PM10的RMSE降低得更多是由于疊加了PM2.5的效果。

        a:地面PM2.5觀測值;b:PM2.5背景場;c:PM2.5同化分析增量;d:PM2.5分析場。

        a:地面PM10觀測值;b:PM10背景場;c:PM10同化分析增量;d:PM10分析場。

        圖5 D02區(qū)域PM2.5和PM10初始場效果檢驗

        2.2 模擬預(yù)報結(jié)果比對分析

        圖6為同化前后PM2.5和PM10的24 h預(yù)報平均濃度值。從圖中可以看到,PM2.5和PM10的實際觀測平均濃度值比Control試驗在初始時刻分別高出約45 μg/m3和70 μg/m3,Da試驗的PM2.5和PM10起始濃度值則分別達到了128 μg/m3和195 μg/m3,同化后兩者的平均濃度更接近觀測值,在模式預(yù)報值普遍低于觀測值的情況下,同化對于初始場的增進效果十分明顯。24 h之內(nèi),同化后的濃度平均值也始終高于模式預(yù)報值,說明同化在整個預(yù)報時段持續(xù)起改進作用;隨著預(yù)報時次的推移,試驗的濃度平均值也隨著模式積分逐漸衰減,PM10的衰減程度較PM2.5更為明顯;后4 h之內(nèi),兩者的濃度值基本達到一致??傮w而言,模式預(yù)報的PM2.5和PM10濃度值基本與觀測曲線走向一致,但預(yù)報值較觀測值整體偏低,這主要取決于模式排放源文件對于模擬結(jié)果的影響。

        圖6 控制試驗和同化試驗PM2.5和PM10的24 h預(yù)報濃度平均值統(tǒng)計

        圖7為同化前后PM2.5和PM10與站點觀測值的RMSE統(tǒng)計結(jié)果,可以看到在起報時刻同化對于PM2.5和PM10均有十分明顯的改進,PM2.5和PM10的RMSE分別平均降低了17.9 μg/m3和24.0 μg/m3。PM2.5的改善程度整體優(yōu)于PM10,改善效果均持續(xù)到了預(yù)報結(jié)束的第24小時,在前18 h改善效果較為明顯,18 h后兩者的RMSE差別逐漸減小,最終趨于相同。圖8為同化前后模式PM2.5和PM10與觀測值的CORR演變曲線,兩者的相關(guān)系數(shù)分別平均提高0.2和0.18。與RMSE演變曲線類似,Da試驗的改進在前20 h比較顯著,其后與Control試驗趨于相同。值得注意的是,無論對于PM2.5還是PM10,Da試驗和Control試驗的相關(guān)系數(shù)在18 h后均降低到0.2以下,而在第10小時左右有一個明顯的峰值,這可能與排放源和大氣環(huán)流的日變化相關(guān)。

        3 結(jié)論

        本文基于WRF-Chem模式和三維變分同化系統(tǒng),利用“全國城市空氣質(zhì)量實時發(fā)布平臺”發(fā)布的實時氣溶膠觀測資料,對華北地區(qū)一次大氣污染過程進行PM2.5和PM10的同化試驗和模擬預(yù)報,并利用統(tǒng)計指標進行檢驗比對分析。結(jié)果表明,資料同化能有效改進背景場模擬偏低的問題,區(qū)域污染分布和強度均發(fā)生改變,與實況更為接近;同化后的初始場,相對于模式背景場,PM2.5和PM10的RMSE平均降低了17.9 μg/m3和24.0 μg/m3,CORR分別平均提高了0.2和0.18。同化對于預(yù)報的改進效果可以持續(xù)到24 h,在24 h內(nèi)氣溶膠平均濃度始終高于模式預(yù)報,但改進效果在前18小時內(nèi)較為突出,隨后效果隨時間迅速降低,其中PM10的衰減程度比PM2.5更為突出,這可能是由于PM10沉降較快,在大氣中停留時間短,故初始場改進對預(yù)報的持續(xù)影響時間也短。

        圖7 控制試驗和同化試驗PM2.5/PM10預(yù)報結(jié)果的RMSE統(tǒng)計

        圖8 控制試驗和同化試驗PM2.5/PM10預(yù)報結(jié)果的CORR統(tǒng)計

        [1] Seinfeld J H, Pandis S N. Atmospheric chemistry and physics: from air pollution to climate change (Second Edition)[M]. Wiley, 1986.

        [2] Kaufman Y J, Tanré D, Boucher O. A satellite view of aerosols in the climate system[J]. Nature, 2002, 419: 215-223.

        [3] 薛文博,王金南,楊金田,等.國內(nèi)外空氣質(zhì)量模型研究進展[J].環(huán)境與可持續(xù)發(fā)展,2013(3):14-20.

        [4] 官元紅,周廣慶,陸維松,等.資料同化方法的理論發(fā)展及應(yīng)用綜述[J].氣象與減災(zāi)研究,2007,30(4):1-8.

        [5] Elbern H, Schmidt H, Ebel A. Variational data assimilation for tropospheric chemistry modeling[J]. Journal of Geophysical Research, 1997, 102: 15967-15985.

        [6] 崔應(yīng)杰,王自發(fā),朱江,等.空氣質(zhì)量數(shù)值模式預(yù)報中資料同化的初步研究[J].氣候與環(huán)境研究,2006(5):616-626.

        [7] Denby B, Schaap M, Segers A, et al. Comparison of data assimilation methods for assessing PM10exceedances on the European scale[J]. Atomspheric Environment, 2008, 42(30): 7122-7134.

        [8] Jiang Z,Liu Z,Wang T,et al. Probing into the impact of 3DVAR assimilation of surface PM10observations over China using process analysis [J]. Journal of Geophysical Research: Atmospheres, 2013, 118(12): 6738-6749.

        [9] Li Z, Zang Z, Li Q B, et al. A three-dimensional variational data assimilation system for multiple aerosol species with WRF/Chem and an application to PM2.5prediction[J]. Atmospheric Chemistry & Physics & Discussions, 2013, 12(8): 4265-4278.

        [10] 韓素芹,馮銀廠,邊海,等.天津大氣污染物日變化特征的WRF-Chem數(shù)值模擬[J].中國環(huán)境科學,2008(9):828-832.

        [11] 章國材.美國WRF模式的進展和應(yīng)用前景[J].氣象,2004,30(12):27-31.

        [12] 周廣強,耿福海,許建明,等.上海地區(qū)臭氧數(shù)值預(yù)報[J].中國環(huán)境科學,2015,36(6):1601-1609.

        [13] Grell G A, Peckham S E, Schmitz R, et al. Fully coupled “online” chemistry within the WRF model [J]. Atmospheric Environment, 2005, 39(37): 6957-6975.

        [14] Ackermann I J, Hass H, Memmesheimer M, et al. Modal aerosol dynamics model for Europe: development and first applications [J]. Atmospheric Environment, 1998, 32(17): 2981-2999.

        [15] Schell B, Ackermann I J, Hass H, et al. Modeling the formation of secondary organic aerosol within a comprehensive air quality model system [J]. Journal of Geophysical Research Atmospheres, 2001, 106(D22): 28275 -28293.

        [16] Zavri R A, Eater R C, Fast J D, et al. Model for simulating aerosol interactions and chemistry(MOSAIC)[J]. Journal of Geophysical Research: Atmospheres(1984~2012), 2008, 113 (D13): 204-218.

        DataAssimilationandSimulationTestsforPM2.5/PM10inAnAtmosphericPollutionProcessoverNorthChina

        HUANG Ran1,2, PAN Xiao-bin1, LI Yi1, HU Fu-jia3

        (1. College of Meteorology and Oceanography, National University of Defense Technology, Nanjing 211101, China; 2. Unit No. 61867, People’s Liberation Army, Beijing 101104, China; 3. Unit No. 95809, People's Liberation Army, Cangzhou 061736, China)

        The PM2.5and PM10aerosol assimilation system was built by using the MOSAIC aerosol variables of WRF-Chem model, and this system was used to carry out the data assimilation and simulation tests for PM2.5/PM10in the atmospheric pollution process over North China on December 29th, 2015. The statistical results indicated that the assimilation could effectively improve the spatial distribution of PM2.5and PM10in initial field, the root-mean-square error of PM2.5and PM10was decreased by 43.7 μg/m3and 71 μg/m3respectively, their correlation coefficient was increased both by 0.52, and the initial location and pollution intensity of heavily polluted region were closer to the actual situation. In addition, the assimilation also could obviously improve the subsequent aerosol forecast, the root-mean-square error of simulated PM2.5and PM10within 24 h was reduced by an average of 17.9 μg/m3and 24.0 μg/m3respectively, and their correlation coefficient was increased by an average of 0.20 and 0.18 separately. However, due to the continuous forcing effect of the emission sources, the improvement effect of initial field assimilation gradually weakened after 18 h.

        WRF-Chem model; Three-dimensional variation; Aerosol; Data assimilation

        2017-06-28

        國家自然科學基金項目(42175128)。

        黃然(1989─),男,重慶人,碩士研究生,研究方向:氣溶膠資料同化與數(shù)值模擬。

        X513

        A

        1001-8581(2017)11-0081-06

        (責任編輯:許晶晶)

        猜你喜歡
        氣溶膠增量觀測
        觀測到恒星死亡瞬間
        軍事文摘(2023年18期)2023-11-03 09:45:42
        提質(zhì)和增量之間的“辯證”
        當代陜西(2022年6期)2022-04-19 12:12:22
        氣溶膠傳播之謎
        “價增量減”型應(yīng)用題點撥
        氣溶膠中210Po測定的不確定度評定
        天測與測地VLBI 測地站周圍地形觀測遮掩的討論
        四川盆地秋季氣溶膠與云的相關(guān)分析
        可觀測宇宙
        太空探索(2016年7期)2016-07-10 12:10:15
        基于均衡增量近鄰查詢的位置隱私保護方法
        電信科學(2016年9期)2016-06-15 20:27:25
        高分辨率對地觀測系統(tǒng)
        太空探索(2015年8期)2015-07-18 11:04:44
        国产在线拍偷自拍偷精品| 中文字幕乱码无码人妻系列蜜桃 | 国内精品91久久久久| 在线视频一区二区观看| 国产三级精品三级在线专区| 午夜精品久久久久久久无码| 蜜桃视频无码区在线观看 | 亚洲成a人v欧美综合天堂| 亚洲日韩一区二区一无码 | 中日韩欧美高清在线播放| 女同国产日韩精品在线| 澳门蜜桃av成人av| 国产亚州精品女人久久久久久| 国产成人av性色在线影院色戒| 97福利视频| 亚洲国产色图在线视频| 中文字幕乱码在线婷婷| 日本免费一区二区三区在线播放| 亚洲一区自拍高清亚洲精品| 粗大的内捧猛烈进出在线视频| 亚洲AV秘 无码一区二区在线| 日本久久精品福利视频| 中文字幕无线码一区二区| 伊人激情av一区二区三区| 国产亚洲欧美在线| 国产在线高清无码不卡| 亚洲精品综合久久中文字幕| 亚洲一区二区三区国产| 国自产拍偷拍精品啪啪一区二区| 欧美粗大无套gay| 天天躁日日操狠狠操欧美老妇| 日本精品啪啪一区二区| 精品露脸熟女区一粉嫩av| 中文字幕乱偷无码av先锋蜜桃| 免费一级肉体全黄毛片| 日本肥老熟妇在线观看| 日本免费大片一区二区三区| 久久精品网站免费观看| 久久久久亚洲av无码专区首jn| 成人免费ā片在线观看| 成人无码无遮挡很H在线播放|