• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    珠江三角洲SO2初始場同化試驗(yàn)研究

    2017-05-23 11:04:18陳懿昂鄧雪嬌高郁東南京信息工程大學(xué)氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室氣候與環(huán)境變化國際合作聯(lián)合實(shí)驗(yàn)室氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心中國氣象局氣溶膠與云降水重點(diǎn)開放實(shí)驗(yàn)室江蘇南京10044中國氣象局廣州熱帶海洋氣象研究所廣東省區(qū)域數(shù)值天氣預(yù)報(bào)重點(diǎn)實(shí)驗(yàn)室廣東廣州510080
    中國環(huán)境科學(xué) 2017年5期
    關(guān)鍵詞:站點(diǎn)尺度背景

    陳懿昂,鄧雪嬌,朱 彬,鄧 濤,高郁東(1.南京信息工程大學(xué)氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室,氣候與環(huán)境變化國際合作聯(lián)合實(shí)驗(yàn)室,氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心,中國氣象局氣溶膠與云降水重點(diǎn)開放實(shí)驗(yàn)室,江蘇 南京 10044;.中國氣象局廣州熱帶海洋氣象研究所,廣東省區(qū)域數(shù)值天氣預(yù)報(bào)重點(diǎn)實(shí)驗(yàn)室,廣東 廣州510080)

    珠江三角洲SO2初始場同化試驗(yàn)研究

    陳懿昂1,2,鄧雪嬌2*,朱 彬1*,鄧 濤2,高郁東2(1.南京信息工程大學(xué)氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室,氣候與環(huán)境變化國際合作聯(lián)合實(shí)驗(yàn)室,氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心,中國氣象局氣溶膠與云降水重點(diǎn)開放實(shí)驗(yàn)室,江蘇 南京 210044;2.中國氣象局廣州熱帶海洋氣象研究所,廣東省區(qū)域數(shù)值天氣預(yù)報(bào)重點(diǎn)實(shí)驗(yàn)室,廣東 廣州510080)

    基于WRF-CMAQ空氣質(zhì)量模式系統(tǒng),采用最優(yōu)插值法(OI)和集合平方根濾波法(EnSRF)對(duì)2013年12月廣東珠江三角洲地區(qū)污染物SO2進(jìn)行初始場同化試驗(yàn).背景場誤差水平分布高值區(qū)主要位于江門一帶,垂直分布在邊界層內(nèi)較大,400m以下基本不變,400m以上隨高度增加而遞減.對(duì)比OI與EnSRF兩種方法同化前后SO2濃度場變化,表明同化觀測資料調(diào)整了污染物濃度場的分布型態(tài),與觀測場更為吻合,兩種方法都能為模式提供與實(shí)際更加接近的初始場.敏感性試驗(yàn)表明最優(yōu)插值法的最優(yōu)水平尺度為 20km;同化站點(diǎn)和檢驗(yàn)站點(diǎn)的均方根誤差平均下降率分別達(dá)到73%和39%.

    資料同化;SO2;最優(yōu)插值法(OI);集合平方根濾波法(EnSRF);珠江三角洲

    近年來,基于CMAQ模式的區(qū)域尺度空氣質(zhì)量預(yù)報(bào)業(yè)務(wù)已在珠江三角洲廣泛地開展,能對(duì)空氣質(zhì)量做出較好預(yù)報(bào)[1-3],但由于模式存在多種不確定性,預(yù)報(bào)的準(zhǔn)確度和時(shí)效性受到不同程度影響.改善模擬效果的途徑之一就是應(yīng)用資料同化方法.

    資料同化是在20世紀(jì)60年代隨著氣象數(shù)值模式發(fā)展而發(fā)展起來,它是通過最優(yōu)估計(jì)算法將觀測與數(shù)值模式信息相結(jié)合以獲得一個(gè)更為真實(shí)的系統(tǒng)狀態(tài),在氣象、海洋領(lǐng)域不斷發(fā)展和完善并已廣泛應(yīng)用.而大氣化學(xué)領(lǐng)域的資料同化研究在20世紀(jì)90年代才開始,從簡單插值方法逐漸發(fā)展到考慮誤差協(xié)方差演變的高級(jí)方法,學(xué)者們采用了不同方法對(duì)模式初始場進(jìn)行優(yōu)化,研究的污染物主要以 O3和氣溶膠顆粒物為主.Lee[4]和 Tombetted等[5]分別采用最優(yōu)插值法(OI)對(duì)韓國和歐洲的 PM10進(jìn)行同化優(yōu)化試驗(yàn),Li等[6]和Jiang等[7]對(duì)WRF-Chem模式中不同氣溶膠方案進(jìn)行了三維變分同化試驗(yàn),Pagowski等[8]對(duì)比了三維變分(3D-Var)和集合卡爾曼濾波(EnKF)對(duì)WRF-Chem模擬O3和PM2.5的優(yōu)化效果.Wu等[9]對(duì)比4種不同方法(OI、EnKF、RRSQRT、4D-Var)對(duì)O3初始場及預(yù)報(bào)的優(yōu)化效果,結(jié)果表明OI同化優(yōu)化初始場的效果最好,而EnKF對(duì)模擬預(yù)報(bào)的優(yōu)化最優(yōu).資料同化的應(yīng)用也從初始條件延伸到了排放源、物理方案參數(shù)化等多種不確定性的修正約束.Elbern等[10]擴(kuò)展EURAD模式的四維變分(4D-Var)同化系統(tǒng)以同時(shí)調(diào)整初始場及物種的排放通量,效果優(yōu)于單純調(diào)整初始場, Vira等[11-12]對(duì) SILAM 化學(xué)傳輸模式進(jìn)行了考慮初始場和首要污染物排放源的同化試驗(yàn).目前大氣化學(xué)資料同化研究是國際上的一個(gè)研究熱點(diǎn)[13-14],中國很多學(xué)者也逐漸加入這個(gè)研究行列中[15-18],對(duì)不同地區(qū)的污染物進(jìn)行了同化分析試驗(yàn).唐曉等[19-20]利用 EnKF對(duì) O3、NOx和VOCs的進(jìn)行初始場及排放源的敏感性實(shí)驗(yàn),探究其對(duì)北京及周邊地區(qū) O3預(yù)報(bào)的改進(jìn)效果,靳璐濱等[21]利用 3D-Var對(duì)青奧會(huì)期間南京地區(qū)PM2.5和PM10進(jìn)行同時(shí)同化試驗(yàn).張金譜等[22]探討集合最優(yōu)插值方法(EnOI)對(duì)珠江三角洲地區(qū)PM10、NO2和SO2初始場的修正能力.

    OI法簡單、計(jì)算量小,但其假設(shè)背景場誤差協(xié)方差是定常.對(duì)于集合平方根濾波法 EnSRF而言,背景場協(xié)方差是隨時(shí)間而變,通過集合預(yù)報(bào)來統(tǒng)計(jì)得出,但其計(jì)算量大.本文嘗試采用 OI法和EnSRF法開展珠江三角洲SO2初始場同化試驗(yàn)研究,結(jié)合地面觀測站 SO2數(shù)據(jù)改進(jìn)CMAQ的初始場,并對(duì)同化站點(diǎn)數(shù)與站點(diǎn)影響范圍的最優(yōu)水平尺度進(jìn)行敏感性分析,以期為珠江三角洲空氣質(zhì)量模式提供更精確的初始場建立同化技術(shù)方法.

    1 所用模式及資料

    本研究所采用的是 WRF-CMAQ[23-24]空氣質(zhì)量模型系統(tǒng),模式采用三重嵌套網(wǎng)格,以廣州為模擬區(qū)域中心,第一區(qū)域?yàn)闁|亞地區(qū),網(wǎng)格數(shù)為182×138,第二區(qū)域包括整個(gè)華南地區(qū),網(wǎng)格數(shù)為98×74,第三區(qū)域?yàn)橹榻侵薜貐^(qū),網(wǎng)格數(shù)為152×110,水平分辨率分別是27、9和3km,同化區(qū)域是第三區(qū)域.排放源用的是MEIC2010年的排放清單.

    圖1 珠江三角洲區(qū)域觀測站點(diǎn)分布Fig.1 Distribution of the Pearl River Delta regional observation stations

    觀測資料來源于廣東省環(huán)境監(jiān)測中心,包括珠江三角洲地區(qū)共58個(gè)觀測站點(diǎn)2013年12月污染物SO2每日逐時(shí)濃度值數(shù)據(jù).其中包括16個(gè)為粵港澳珠江三角洲區(qū)域空氣監(jiān)測網(wǎng)絡(luò)站點(diǎn),分別為廣州的天湖(TH)、磨碟沙(MDS)、竹洞(ZD)、麓湖(LH)、深圳荔園(LY)、珠海唐家(TJ)、佛山的惠景城(HJC)、金桔咀(JJJ)、江門的端芬(DF)、東湖(DH)、肇慶城中子站(CZZ)、惠州的下埔(XP)、金果灣(JGW)、西角(XJ)、東莞南城元嶺(NCYL)、中山紫馬嶺(ZML),站點(diǎn)位置如圖1所示.對(duì)于同化站點(diǎn)的選擇,要綜合考慮站點(diǎn)覆蓋范圍,數(shù)據(jù)代表性和數(shù)據(jù)質(zhì)量等因素,這里以粵港澳珠江三角洲區(qū)域空氣監(jiān)測網(wǎng)絡(luò)站點(diǎn)為主,再增加部分站點(diǎn)作為同化站點(diǎn).參考張金譜的研究[22],試驗(yàn)中設(shè)置較多檢驗(yàn)站點(diǎn)以避免代表性不足對(duì)同化整體效果的錯(cuò)誤分析.

    2 方法介紹

    本文選用OI法和EnSRF法,兩者最主要的區(qū)別在于背景場誤差的計(jì)算.OI法是由Eliassen[25]提出的,在統(tǒng)計(jì)意義上通過使分析方差最小來得到最優(yōu)分析值,其計(jì)算復(fù)雜度和計(jì)算量較小,背景場誤差定常.EnSRF法是由Whitaker和 Hamill[26]提出的不加入觀測擾動(dòng)的集合卡爾曼濾波,通過有限樣本來計(jì)算背景場誤差協(xié)方差,計(jì)算量較大,背景場誤差隨時(shí)間而變.

    兩種方法的基本公式如下:

    EnSRF式中:

    xa為同化后分析場;xb為背景場;xf12和 xf24分別表示12h和24h預(yù)報(bào)場;yo為觀測場;H為觀測算子,將狀態(tài)量從模式空間轉(zhuǎn)換到觀測空間;Bb表示背景場誤差協(xié)方差矩陣;R表示觀測場誤差協(xié)方差矩陣,這里不考慮站點(diǎn)間的相關(guān)性,所以R為對(duì)角矩陣,以觀測值的 10%作為觀測誤差[17];K 為結(jié)合了背景誤差協(xié)方差和觀測誤差協(xié)方差得到的權(quán)重算子;N為樣本個(gè)數(shù),本文中樣本數(shù)為30個(gè);上標(biāo)“T”和“-1”分別表示矩陣的轉(zhuǎn)置和求逆,上標(biāo)“-”和“’”分別表示樣本平均和樣本偏差.

    背景誤差協(xié)方差是資料同化的核心部分,它決定著觀測信息在同化分析中的傳播,對(duì)于最優(yōu)插值法,這里采用 National Meteorological Center(NMC)[27]的方法來計(jì)算背景場誤差協(xié)方差,即利用同一時(shí)刻的兩個(gè)不同預(yù)報(bào)時(shí)效的預(yù)報(bào)結(jié)果之差來表示背景場誤差,模式資料包括了2013年12月01~30日12h和24h的0:00預(yù)報(bào)場,將這兩個(gè)預(yù)報(bào)場作差以得到背景誤差場.將其分解為誤差標(biāo)準(zhǔn)差矩陣和相關(guān)系數(shù)矩陣兩部分,對(duì)于相關(guān)系數(shù)矩陣,又可分水平和垂直相關(guān),這里進(jìn)行地面同化,所以僅考慮水平相關(guān)系數(shù) Cx,采用一維的高斯分布函數(shù)來表示,假設(shè)水平各向同性,其中 xΔ 表示兩點(diǎn)間的水平距離,L表示相關(guān)系數(shù)的消減尺度.而對(duì)于集合平方根濾波法,則是利用集合樣本來計(jì)算,樣本的生成方法參考Evensen[28]提出的三維隨機(jī)場生成方法,生成均值為0,方差為1的隨機(jī)擾動(dòng)場,將初始值加上一定大小的擾動(dòng)形成N個(gè)初始樣本,進(jìn)行一定時(shí)間的集合預(yù)報(bào),得到背景場.

    3 結(jié)果與討論

    3.1 模式性能評(píng)估

    先對(duì)WRF-CMAQ在所選時(shí)間段的模擬情況進(jìn)行效果評(píng)估,對(duì)2013年12月1~25日SO2的日均濃度在全部站點(diǎn)的模擬效果統(tǒng)計(jì)分析,從平均標(biāo)準(zhǔn)誤差(NMB)上看,其極值體現(xiàn)了模式在不同站點(diǎn)模擬偏離程度的極端情況,可知同時(shí)存在偏高和偏低的情況,均值為 61%,說明模式對(duì)SO2的模擬整體偏高.相關(guān)系數(shù)(R)均值達(dá)到0.75,但仍有站點(diǎn)僅為 0.34,說明模式對(duì)不同區(qū)域的模擬具有不穩(wěn)定性.按照不同區(qū)域?qū)φ军c(diǎn)進(jìn)行評(píng)估,也可以發(fā)現(xiàn)模式對(duì)于廣佛地區(qū)的模擬較好,對(duì)于珠海深圳等地的模擬效果較差.

    表1 SO2日均濃度模擬效果統(tǒng)計(jì)分析Table 1 Statistical analyses of modelling performance on daily mean SO2

    3.2 背景場分析

    圖2 SO2濃度在模式第一層分布Fig.2 Distribution of the SO2concentration at the first model level

    圖2a是SO2的第一層月平均濃度分布,可見在佛山、江門和中山3市交界處有一個(gè)高值區(qū),深圳市也有一個(gè)高值區(qū),而通過 NMC方法統(tǒng)計(jì)得到 SO2的背景場誤差標(biāo)準(zhǔn)差的區(qū)域分布(圖 2b)與濃度場分布類似,但是高值區(qū)的對(duì)應(yīng)不是很好,誤差標(biāo)準(zhǔn)差的高值區(qū)主要位于區(qū)域南部江門一帶.OI法的背景場誤差為靜態(tài),不隨時(shí)間變化,而EnSRF法的背景場誤差是動(dòng)態(tài)的,隨模式積分而更新,理論上要優(yōu)于 OI法的背景場誤差.圖3是通過集合計(jì)算出來的12月9~11日的誤差標(biāo)準(zhǔn)差,可見其分布型態(tài)是隨時(shí)間變化的.背景場整體濃度大小10日>11日>9日,整體的誤差標(biāo)準(zhǔn)差大小也是10日>11日>9日,即誤差變化趨勢與濃度變化一致.而兩種方法計(jì)算所得的誤差高值區(qū)都主要位于區(qū)域南部江門一帶,說明模式對(duì)此處的模擬不確定性較大.

    圖3 12月9~11日的誤差標(biāo)準(zhǔn)差分布Fig.3 Distribution of the standard error from Dec 9 to Dec 11

    圖 4是背景場誤差標(biāo)準(zhǔn)差和平均濃度的垂直廓線,從中可見SO2的平均濃度在400m以下基本維持不變,從 400m往上,濃度值隨高度的升高迅速遞減,到4km左右的高度濃度值幾乎為0,說明SO2在垂直方向400m內(nèi)混合比較均勻,分布范圍大概能達(dá)到4km左右的高度,這與其物理化學(xué)性質(zhì)的穩(wěn)定性有關(guān),是影響它傳播尺度的因素之一.而背景場誤差標(biāo)準(zhǔn)差也有相似的變化趨勢,誤差標(biāo)準(zhǔn)差越大,說明模式對(duì)SO2模擬的不確定性越大,可以看出誤差標(biāo)準(zhǔn)差在邊界層內(nèi),特別在400m以下比較大,這是由于SO2主要來自近地面的各種源排放,邊界層內(nèi)湍流運(yùn)動(dòng)強(qiáng)且復(fù)雜,而且 SO2的化學(xué)性質(zhì)活潑,易與其它物質(zhì)發(fā)生化學(xué)反應(yīng)生成二次污染物,這些物理化學(xué)變化都影響著 SO2濃度的變化,所以其在邊界層內(nèi)的模擬不確定性較大.

    圖4 SO2垂直廓線圖Fig.4 Vertical profile of SO2

    3.3 OI與EnSRF同化效果分析

    圖5 不同方法同化前后檢驗(yàn)站點(diǎn)的MAGE(上)與RMSE(下)對(duì)比Fig.5 Comparisons of MAGE(top) and RMSE(bottom) before and after assimilation using different methods

    采用OI法和EnSRF法分別對(duì)12月9~11日 0:00模式初始場進(jìn)行同化試驗(yàn),共同化了 31個(gè)觀測站點(diǎn)(圖1中●+○)數(shù)據(jù).從統(tǒng)計(jì)上看,2種方法同化后站點(diǎn)誤差均有所下降,平均絕對(duì)偏差(MAGE)和均方根誤差(RMSE)平均分別從33和39μg/m3下降到17和21μg/m3左右(表2),EnSRF法的整體效果要略好于OI法,但在不同區(qū)域站點(diǎn)的同化效果各有優(yōu)劣.圖5是2種方法同化前后站點(diǎn)的MAGE和RMSE對(duì)比,其中2種方法對(duì)江門地區(qū)的同化效果都比較差,這與江門地區(qū)所選同化站點(diǎn)區(qū)域代表性有關(guān).OI法對(duì)佛山地區(qū)站點(diǎn)的同化效果要優(yōu)于 EnSRF法,OI法同化后在廣州地區(qū)部分站點(diǎn)出現(xiàn) RMSE反而增大的現(xiàn)象,EnSRF則不會(huì)出現(xiàn)這種情況,其余大部分地區(qū)EnSRF法優(yōu)化效果比 OI法略好.從濃度區(qū)域分布圖可見,模式對(duì)SO2濃度模擬整體偏高,背景場從中部往西南部存在高值區(qū),而且從多個(gè)不同時(shí)次的背景場中發(fā)現(xiàn)污染物分布型態(tài)基本一致,高值區(qū)基本也都位于這一區(qū)域,這可能與排放源、氣象場的風(fēng)向風(fēng)速模擬以及模式自身誤差相關(guān).通過 2種不同方法進(jìn)行同化后,濃度場的分布均發(fā)生了變化,原高值區(qū)域的濃度有了不同程度的下降,同化后高值區(qū)主要位于廣州北部及江門與佛山交界處,與觀測場濃度更加貼近.但由于2種方法背景場誤差的大小存在差異,使得同化對(duì)背景場的調(diào)整力度不一樣,如 10日的肇慶區(qū)域(圖6),從背景場誤差上(圖2b與圖3b)EnSRF法的比OI法的大,導(dǎo)致分析增量比較大,使得EnSRF法同化后此區(qū)域濃度變得比較低,從統(tǒng)計(jì)上看在這個(gè)區(qū)域EnSRF法的整體效果要優(yōu)于OI法.上述結(jié)果表明,OI和EnSRF2種方法都能提供一個(gè)與實(shí)際更接近的初始場.

    表2 不同方法同化結(jié)果統(tǒng)計(jì)(μg/m3)Table 2 Statistics of the different assimilation methods (μg/m3)

    圖6 SO2濃度場(a背景場 b OI分析場 c EnSRF分析場μg/m3)Fig.6 SO2concentration fields (a background field b OI analysis field c EnSRF analysis field μg/m3)

    3.4 敏感性試驗(yàn)分析

    上述試驗(yàn)中EnSRF法與OI法優(yōu)化效果相近,且OI法計(jì)算量較小,所以在此采用OI法進(jìn)行同化敏感性試驗(yàn)分析.同化試驗(yàn)中需要考慮同化站點(diǎn)的選擇及水平相關(guān)尺度L的設(shè)置,因?yàn)檎军c(diǎn)分布疏密及區(qū)域代表性好壞和水平尺度大小是否合適都對(duì)同化效果的優(yōu)劣有所影響,所以下面分別對(duì)這兩個(gè)方面進(jìn)行了敏感性試驗(yàn).

    在選擇同化站點(diǎn)時(shí)進(jìn)行了 3次試驗(yàn)(這里水平相關(guān)尺度L均取20km),以粵港澳珠江三角洲洲區(qū)域監(jiān)測網(wǎng)絡(luò)的16個(gè)站點(diǎn)為主,考慮站點(diǎn)的覆蓋范圍,分別同化了16個(gè)(圖1中●)、31個(gè)(圖1中●+○)和全部58個(gè)(圖1中●+○+△)站點(diǎn)的觀測.以3次試驗(yàn)中相同的同化站點(diǎn)(即圖1中●)和檢驗(yàn)站點(diǎn)(即圖1中△)進(jìn)行統(tǒng)計(jì)對(duì)比(全部站點(diǎn)參與同化時(shí)沒有檢驗(yàn)站點(diǎn),圖7b中圖例58的數(shù)據(jù)可用來對(duì)比站點(diǎn)同化與未同化的差別),由圖7a可見,3次試驗(yàn)中同化站點(diǎn)的均方根誤差(RMSE)均比未同化時(shí)有較大減小.從統(tǒng)計(jì)上看,隨著站點(diǎn)數(shù)的增加,同化站點(diǎn)的RMSE越大,但這不表明同化站點(diǎn)數(shù)越多所得到的分析場反而不好.例如同化站點(diǎn)相鄰比較近的時(shí)候,同化更多站點(diǎn)時(shí)所得分析場,是受更多觀測值影響的加權(quán)平均,所以統(tǒng)計(jì)所得的RMSE會(huì)比同化少數(shù)站點(diǎn)時(shí)大.增加觀測資料,如31個(gè)站點(diǎn)的試驗(yàn),增加了站點(diǎn)覆蓋面,有利于站點(diǎn)附近的改進(jìn),而對(duì)于58個(gè)站點(diǎn)的試驗(yàn),站點(diǎn)分布較密區(qū)域雖然受到共同影響,RMSE有所增加,但同時(shí)也可以減小觀測資料代表性不好的影響.同化站點(diǎn)選擇的最優(yōu)化,除了考慮覆蓋面和區(qū)域代表性,還需通過預(yù)報(bào)場進(jìn)一步檢驗(yàn),這也是下一步要開展的研究.不同同化站點(diǎn)間會(huì)相互影響,特別當(dāng)站點(diǎn)數(shù)多且站點(diǎn)間距離較近時(shí)影響更深,所以需通過設(shè)置合適水平相關(guān)尺度,使得站點(diǎn)間的影響合理化,從而使同化效果最優(yōu)化.

    水平相關(guān)尺度L的大小決定著同化站點(diǎn)的影響范圍,如果設(shè)置過小,觀測分布稀疏地區(qū)的很多區(qū)域不能得到有效更新,不能充分利用觀測數(shù)據(jù),如果設(shè)置過大,則觀測站影響區(qū)域偏大,可能生成與事實(shí)不符的分析場,反而使得同化效果下降.本文選取了不同水平相關(guān)尺度對(duì)2013年12月每日0:00(世界時(shí))的SO2預(yù)報(bào)場進(jìn)行同化試驗(yàn),同化站點(diǎn)與上面同化站點(diǎn)數(shù)為31的站點(diǎn)一致,再對(duì)同化站點(diǎn)和檢驗(yàn)站點(diǎn)進(jìn)行統(tǒng)計(jì)對(duì)比,以選出此試驗(yàn)時(shí)段的最優(yōu)水平尺度,因?yàn)椴煌军c(diǎn)模式值與觀測值的偏差有正有負(fù),所以這里采用平均絕對(duì)偏差(MAGE)和均方根誤差(RMSE)進(jìn)行統(tǒng)計(jì).

    圖7 同化不同數(shù)目觀測值的RMSE對(duì)比Fig.7 Comparisons of RMSE with assimilating different number of observations

    表3 不同水平相關(guān)尺度同化結(jié)果統(tǒng)計(jì)(μg/m3)Table 3 Statistics of the assimilation result with different horizontal correlation scale (μg/m3)

    從表 3中可見,采用不同尺度得到的分析場的MAGE和RMSE均小于背景場.對(duì)于同化站點(diǎn),隨著相關(guān)尺度的增加,MAGE和 RMSE均變大,這可能是因?yàn)殡S相關(guān)尺度的增大,觀測站點(diǎn)能影響的范圍變大,加深了不同站點(diǎn)間的相互影響,使得同化站點(diǎn)受到更多站點(diǎn)的影響,當(dāng)站點(diǎn)間差異較大時(shí),就使得同化效果下降.對(duì)于檢驗(yàn)站點(diǎn),背景場的RMSE約為29.5μg/m3,相關(guān)尺度為10km、20km和30km時(shí)的同化效果差別不大,RMSE減少到17.7μg/m3左右,其中20km的RMSE較小,為17.6μg/m3,偏差也隨著相關(guān)尺度的變大而增加.若水平尺度取 10km,則同化站點(diǎn)調(diào)整的區(qū)域范圍較小,如江門市觀測站點(diǎn)分布稀疏的這種情況,大部分區(qū)域不能得到有效的更新,所以綜合上面結(jié)果取這次試驗(yàn)的最優(yōu)水平尺度為20km.

    表4 同化站點(diǎn)及檢驗(yàn)站點(diǎn)同化前后RMSE對(duì)比Table 4 Comparisons of RMSE before and after assimilation at assimilation sites and validation sites

    圖8 同化站點(diǎn)同化前后RMSE對(duì)比Fig.8 Comparisons of RMSE before and after assimilation at assimilation sites

    以水平尺度為20km且同化了31個(gè)站點(diǎn)的同化結(jié)果來進(jìn)行分析.圖8展示了SO2在各個(gè)同化站點(diǎn)同化前后的RMSE對(duì)比情況,可見RMSE均有大幅度的下降,最低也有31%的降幅,其中端芬站的RMSE下降比例最高,從背景場誤差分布上可以看出,端芬站所在江門區(qū)域誤差較大,且觀測誤差較小,所以同化過程對(duì)觀測值比較信任,而下降比例最低的佛山三水監(jiān)測站剛好相反,位于源排放較強(qiáng)的區(qū)域,試驗(yàn)時(shí)段的觀測值較大,按觀測值的10%為觀測誤差的設(shè)置,觀測誤差較大,而背景誤差較小,導(dǎo)致同化效果較差,這也表明正確估算背景誤差和觀測誤差的重要性.

    圖 9是 SO2在各個(gè)檢驗(yàn)站點(diǎn)同化前后的RMSE對(duì)比情況,可見檢驗(yàn)站點(diǎn)的優(yōu)化效果小于同化站點(diǎn),其中大部分站點(diǎn)具有不同程度的下降,但有個(gè)別站點(diǎn)出現(xiàn)了同化后 RMSE反而上升的現(xiàn)象,以江門地區(qū)站點(diǎn)尤為明顯,對(duì)比江門地區(qū)同化的東湖站與其周圍檢驗(yàn)站點(diǎn)觀測均值可以發(fā)現(xiàn),周圍檢驗(yàn)站點(diǎn)濃度偏高,而東湖站的濃度值較低,且它的同化效果較好,所以將背景場這個(gè)區(qū)域濃度值調(diào)整得偏低,導(dǎo)致了其他站點(diǎn)的模式值與觀測值的誤差變大,這也說明在這段時(shí)間的東湖站區(qū)域代表性較差.圖 10是同化前后同化站點(diǎn)SO2模式值與觀測值的散點(diǎn)圖,站點(diǎn)在同化后相關(guān)性有較大提高.圖 10b右下角的離散點(diǎn)位于佛山三水地區(qū),正如前面分析,佛山區(qū)域背景場誤差估算較小,而觀測值較大使得觀測誤差較大,所以分析對(duì)模式結(jié)果較信任,使得同化觀測雖有所優(yōu)化,但效果較弱.

    圖9 檢驗(yàn)站點(diǎn)同化前后RMSE對(duì)比Fig.9 Comparisons of RMSE before and after assimilation at validation sites

    本研究還存在一些不足之處,需要對(duì)觀測誤差進(jìn)行更精確的估算,這里只是對(duì)EnSRF方法的初步嘗試,擾動(dòng)生成樣本方式及很多相關(guān)參數(shù)如集合樣本數(shù),局地化尺度及資料同化的質(zhì)量控制等還需進(jìn)行進(jìn)一步敏感性試驗(yàn)來最優(yōu)化.需要進(jìn)行更長時(shí)間的同化試驗(yàn)來比較兩種同化方法的效果.

    圖10 同化前后模式值與觀測值散點(diǎn)圖Fig.10 Scatter diagram of model value and observation value before and after assimilation

    4 結(jié)論

    4.1 基于珠江三角洲地區(qū)SO21個(gè)月的預(yù)報(bào)結(jié)果和集合預(yù)報(bào)結(jié)果統(tǒng)計(jì)分析表明背景場誤差高值區(qū)主要位于江門一帶,集合預(yù)報(bào)所得背景場誤差的變化趨勢與濃度場變化趨勢一致.模式對(duì)邊界層內(nèi)特別400m以下SO2的模擬不確定性較大,主要跟源排放、氣象條件及SO2自身物理化學(xué)性質(zhì)有關(guān).

    4.2 集合平方根濾波對(duì)SO2初始場的優(yōu)化效果與最優(yōu)插值法的結(jié)果相近,同化觀測資料對(duì)模式濃度值有不同程度的調(diào)整,改變 SO2濃度場的分布,使其與觀測場更為吻合,表明2種方法都具有改善模式初始場的作用.

    4.3 通過敏感性試驗(yàn)確定了此次試驗(yàn)的最優(yōu)水平尺度為20km.SO2初始場經(jīng)OI法同化所得分析值更接近觀測值,站點(diǎn)的均方根誤差(RMSE)和平均絕對(duì)偏差(MAGE)均有所下降,其中同化站點(diǎn)和檢驗(yàn)站點(diǎn)的RMSE平均分別下降73%和39%.從同化不同數(shù)目觀測資料試驗(yàn)的統(tǒng)計(jì)上看,受到更多站點(diǎn)的共同影響,同化站點(diǎn)的RMSE隨著同化站點(diǎn)數(shù)的增多而增大,但相比于未同化時(shí)的RMSE均有所優(yōu)化.

    [1] 鄧雪嬌,鄧 濤,吳 兌,等.珠江三角洲空氣質(zhì)量與能見度數(shù)值預(yù)報(bào)模式系統(tǒng) [J]. 廣東氣象, 2010,32(4):18-22.

    [2] 鄧 濤,吳 兌,鄧雪嬌,等.珠江三角洲空氣質(zhì)量暨光化學(xué)煙霧數(shù)值預(yù)報(bào)系統(tǒng) [J]. 環(huán)境科學(xué)與技術(shù), 2013,36(4):62-68.

    [3] 鄧 濤,鄧雪嬌,吳 兌,等.珠江三角洲灰霾數(shù)值預(yù)報(bào)模式與業(yè)務(wù)運(yùn)行評(píng)估 [J]. 氣象科技進(jìn)展:英文版, 2012,2(6):38-44.

    [4] Lee E H, Ha J C, Lee S S, et al. PM10data assimilation over south Korea to Asian dust forecasting model with the optimal interpolation method [J]. Asia-Pacific Journal of Atmospheric Sciences, 2013,49(1):73-85.

    [5] Tombette M, Mallet V, Sportisse B, et al. PM10data assimilation over Europe with the optimal interpolation method [J]. Atmospheric Chemistry & Physics, 2009,9(1):57-70.

    [6] 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, 2012,12(5):4265-4278.

    [7] 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.

    [8] Pagowski M, Grell G A. Experiments with the assimilation of fine aerosols using an ensemble Kalman filter [J]. Journal of Geophysical Research: Atmospheres (1984–2012), 2012, 117(D21).

    [9] Wu L, Mallet V, Bocquet M, et al. A comparison study of data assimilation algorithms for ozone forecasts [J]. Journal of Geophysical Research: Atmospheres (1984–2012), 2008,113 (D20).

    [10] Elbern H, Strunk A, Schmidt H, et al. Emission rate and chemical state estimation by 4-dimensional variational inversion [J]. Atmospheric Chemistry & Physics Discussions, 2007,7(14): 3749-3769.

    [11] Vira J, Sofiev M. On variational data assimilation for estimating the model initial conditions and emission fluxes for short-term forecasting of SOx concentrations [J]. Atmospheric Environment,2012,46(1):318-328.

    [12] Vira J, Sofiev M. Assimilation of surface NO2and O3observations into the SILAM chemistry transport model [J]. Geoscientific Model Development, 2015,8(2):191-203.

    [13] Bocquet M, Elbern H, Eskes H, et al. Data assimilation in atmospheric chemistry models: current status and future prospects for coupled chemistry meteorology models [J]. Atmospheric Chemistry & Physics Discussions, 2014,15(23):32233-32323.

    [14] Zhang Y, Bocquet M, Mallet V, et al. Real-time air quality forecasting, Part II: State of the science, current research needs, and future prospects [J]. Atmospheric Environment, 2012,60(6): 656–676.

    [15] Niu T, Gong S L, Zhu G F, et al. Data assimilation of dust aerosol observations for the CUACE/dust forecasting system [J]. Atmospheric Chemistry and Physics, 2008,8(13):3473-3482.

    [16] Dai T, Schutgens N A J, Goto D, et al. Improvement of aerosol optical properties modeling over Eastern Asia with MODIS AOD assimilation in a global non-hydrostatic icosahedral aerosol transport model [J]. Environmental Pollution, 2014,195:319-329.

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

    [18] 白曉平,李 紅,方 棟,等.資料同化方法在空氣污染數(shù)值預(yù)報(bào)中的應(yīng)用研究 [J]. 環(huán)境科學(xué), 2008,29(2):283-289.

    [19] Tang X, Zhu J, Wang Z F, et al. Improvement of ozone forecast over Beijing based on ensemble Kalman filter with simultaneous adjustment of initial conditions and emissions [J]. Atmospheric Chemistry & Physics, 2011,11(3):7811-7849.

    [20] 唐 曉,朱 江,王自發(fā),等.基于集合卡爾曼濾波的區(qū)域臭氧資料同化試驗(yàn) [J]. 環(huán)境科學(xué)學(xué)報(bào), 2013,33(3):796-805.

    [21] 靳璐濱,臧增亮,潘曉濱,等.PM2.5和 PM2.5~10資料同化及在南京青奧會(huì)期間的應(yīng)用試驗(yàn) [J]. 中國環(huán)境科學(xué), 2016,36(2):331-341.

    [22] 張金譜,胡嘉鏜,王雪梅.集合最優(yōu)插值同化方法在珠江三角洲空氣質(zhì)量模擬中的初步應(yīng)用 [J]. 環(huán)境科學(xué)學(xué)報(bào), 2014,34(3):558- 566.

    [23] Skamarock W C, Klemp J B, Dudhia J, et al. A description of the Advanced Research WRF version 3 [Z]. NCAR Technical Note:NCAR/TN-475+STR, 2008.

    [24] Byun D W, Ching J K S, Byun D W, et al. Science Algorithms of the EPA Models - 3Community Multiscale Air Quality (CMAQ) Modeling System [J]. Environmental Protection Agency Office of Research & Development, 1999.

    [25] Eliassen A. Provisional Report on Calculation of Spatial Covariance and Autocorrelation of Pressure Field [J]. Space Science Reviews - SPACE SCI REV, 1954:5.

    [26] Whitaker J S, Hamill T M. Ensemble Data Assimilation without Perturbed Observations [J]. Monthly Weather Review, 2002, 130(7):1913-1924.

    [27] Parrish D F, Derber J C. The National Meteorological Center's spectral statistical-interpolation analysis system [J]. Monthly Weather Review, 1992,120(8):1747-1763.

    [28] Evensen G. The Ensemble Kalman Filter: theoretical formulation and practical implementation [J]. Ocean Dynamics, 2003,53(4): 343-367.

    Data assimilation experiment on SO2initial conditions in the Pearl River Delta.

    CHEN Yi-ang1,2, DENG Xue-jiao2*, ZHU Bin1*, DENG Tao2, GAO Yu-dong2(1.Key Laboratory of Meteorological Disaster, Ministry of Education, Joint International Research Laboratory of Climate and Environment Change, Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters, Key Laboratory for Aerosol-Cloud-Precipitation of China Meteorological Administration, Nanjing University of Information Science & Technology, Nanjing 210044, China;2.Guangdong Provincial Key Laboratory of Regional Numerical Weather Prediction, Institute of Tropical and Marine Meteorology, China Meteorological Administration, Guangzhou 510080, China). China Environmental Science, 2017,37(5):1610~1619

    Based on the WRF-CMAQ air quality model, the pollutant SO2in the Pearl River Delta Region in December 2013 was assimilated to optimize the initial conditions using the optimal interpolation approach (OI) and the ensemble square root filter (EnSRF). The high values of the background error were mainly located in Jiangmen region in horizontal direction and were larger within the boundary layer in vertical direction. The background error was nearly constant below 400m and decreased with height above 400m. By comparing the SO2concentration fields using assimilation with those not using assimilation, the results showed that assimilation could adjust the distribution pattern of the pollutant and make it more consistent with the observation field. Both assimilated methods could offer an initial field closer to the true situation. The sensitivity test showed that the optimal horizontal scale of the optimal interpolation method was 20km. The root mean square error decreasing percentage between the assimilation sites and the verification sites reached 73% and 39%, respectively. With the number of the assimilation site increasing, the optimization of the assimilation site had a declining trend.

    data assimilation;SO2;optimal interpolation (OI);ensemble square root filter (EnSRF);Pearl River Delta Region

    X511

    A

    1000-6923(2017)05-1610-10

    陳懿昂(1991-),男,廣東汕頭人,南京信息工程大學(xué)碩士研究生,主要從事大氣化學(xué)大氣環(huán)境方向研究.

    2016-09-18

    國家自然科學(xué)基金資助項(xiàng)目(41475105);國家科技支撐計(jì)劃課題(2014BAC16B06);科技部公益性(氣象)行業(yè)項(xiàng)目(GYHY201306042);中國氣象局氣候變化專項(xiàng)項(xiàng)目(CCSF201531);廣東省科技計(jì)劃項(xiàng)目(2015A020215020);廣東省氣象局科技創(chuàng)新團(tuán)隊(duì)計(jì)劃項(xiàng)目(201506)

    * 責(zé)任作者, 鄧雪嬌, 研究員, dxj@grmc.gov.cn; 朱 彬, 教授, binzhu@nuist.edu.cn

    猜你喜歡
    站點(diǎn)尺度背景
    “新四化”背景下汽車NVH的發(fā)展趨勢
    《論持久戰(zhàn)》的寫作背景
    財(cái)產(chǎn)的五大尺度和五重應(yīng)對(duì)
    基于Web站點(diǎn)的SQL注入分析與防范
    電子制作(2019年14期)2019-08-20 05:43:42
    2017~2018年冬季西北地區(qū)某站點(diǎn)流感流行特征分析
    首屆歐洲自行車共享站點(diǎn)協(xié)商會(huì)召開
    中國自行車(2017年1期)2017-04-16 02:53:52
    晚清外語翻譯人才培養(yǎng)的背景
    怕被人認(rèn)出
    宇宙的尺度
    太空探索(2016年5期)2016-07-12 15:17:55
    9
    侵犯人妻中文字幕一二三四区| 热re99久久国产66热| 色视频在线一区二区三区| 99国产精品一区二区三区| 少妇粗大呻吟视频| 亚洲精品美女久久久久99蜜臀| 亚洲一区二区三区欧美精品| 免费久久久久久久精品成人欧美视频| 夜夜骑夜夜射夜夜干| 午夜精品国产一区二区电影| 一二三四社区在线视频社区8| 丝袜人妻中文字幕| 免费av中文字幕在线| 两性午夜刺激爽爽歪歪视频在线观看 | 十八禁高潮呻吟视频| 国产精品二区激情视频| 国产欧美日韩一区二区精品| 啦啦啦在线免费观看视频4| 亚洲中文字幕日韩| 中文字幕人妻丝袜一区二区| 欧美亚洲 丝袜 人妻 在线| 人妻久久中文字幕网| 我要看黄色一级片免费的| 99久久人妻综合| 成年人免费黄色播放视频| 国产色视频综合| 成人18禁高潮啪啪吃奶动态图| 国产精品1区2区在线观看. | 黑人巨大精品欧美一区二区mp4| 日韩电影二区| 1024视频免费在线观看| 丰满迷人的少妇在线观看| 国产亚洲av高清不卡| 日韩欧美一区视频在线观看| 肉色欧美久久久久久久蜜桃| 久久久国产一区二区| 欧美国产精品va在线观看不卡| 97在线人人人人妻| 91av网站免费观看| 中文字幕色久视频| 国产精品久久久久成人av| 在线观看人妻少妇| 日本a在线网址| 老司机午夜福利在线观看视频 | 久久国产精品影院| 午夜免费观看性视频| 久久精品国产亚洲av高清一级| 国产精品亚洲av一区麻豆| 一区二区三区乱码不卡18| 日本五十路高清| 国产成人精品在线电影| 亚洲成国产人片在线观看| 日韩免费高清中文字幕av| 日本五十路高清| 99精国产麻豆久久婷婷| 熟女少妇亚洲综合色aaa.| 在线十欧美十亚洲十日本专区| av在线老鸭窝| 精品免费久久久久久久清纯 | 亚洲国产av影院在线观看| 91精品国产国语对白视频| 天天躁夜夜躁狠狠躁躁| 日韩欧美免费精品| 爱豆传媒免费全集在线观看| 9色porny在线观看| 老司机影院毛片| 少妇被粗大的猛进出69影院| 亚洲国产看品久久| 99久久精品国产亚洲精品| 国产男人的电影天堂91| 天堂8中文在线网| 青青草视频在线视频观看| 国产在视频线精品| 9191精品国产免费久久| 久久国产亚洲av麻豆专区| 欧美日韩黄片免| 在线观看人妻少妇| 考比视频在线观看| av视频免费观看在线观看| 另类亚洲欧美激情| 亚洲精品国产精品久久久不卡| 精品国产一区二区三区久久久樱花| 国产精品av久久久久免费| 日韩制服骚丝袜av| 2018国产大陆天天弄谢| 亚洲国产欧美在线一区| 色视频在线一区二区三区| 亚洲中文日韩欧美视频| 一本久久精品| 99久久人妻综合| 成人av一区二区三区在线看 | 欧美亚洲 丝袜 人妻 在线| 亚洲av男天堂| 两人在一起打扑克的视频| 婷婷丁香在线五月| 成人三级做爰电影| 中文字幕色久视频| 亚洲精品国产色婷婷电影| 美女福利国产在线| 黄色视频在线播放观看不卡| 纯流量卡能插随身wifi吗| 搡老岳熟女国产| 亚洲男人天堂网一区| 亚洲性夜色夜夜综合| 一级,二级,三级黄色视频| 视频区图区小说| 国产精品欧美亚洲77777| 国产成人精品久久二区二区免费| 91字幕亚洲| 人人妻人人澡人人爽人人夜夜| 亚洲伊人久久精品综合| 一本综合久久免费| 韩国精品一区二区三区| 国产野战对白在线观看| 法律面前人人平等表现在哪些方面 | 色婷婷av一区二区三区视频| 精品少妇一区二区三区视频日本电影| 夫妻午夜视频| 两个人看的免费小视频| 国产在视频线精品| 狂野欧美激情性bbbbbb| 一区二区三区激情视频| 999精品在线视频| 国产又爽黄色视频| av有码第一页| 亚洲精品日韩在线中文字幕| 99久久精品国产亚洲精品| 两性午夜刺激爽爽歪歪视频在线观看 | 丝袜美腿诱惑在线| 中文字幕最新亚洲高清| 女人被躁到高潮嗷嗷叫费观| 人成视频在线观看免费观看| 久久影院123| 纯流量卡能插随身wifi吗| 天天添夜夜摸| 欧美大码av| 12—13女人毛片做爰片一| 国产免费现黄频在线看| 自拍欧美九色日韩亚洲蝌蚪91| 又紧又爽又黄一区二区| 视频区欧美日本亚洲| 爱豆传媒免费全集在线观看| 亚洲国产欧美日韩在线播放| 日韩视频一区二区在线观看| 黑人猛操日本美女一级片| 超碰97精品在线观看| 国产精品一区二区在线观看99| 狂野欧美激情性xxxx| 亚洲精品第二区| 欧美日本中文国产一区发布| 国产精品免费大片| 韩国高清视频一区二区三区| 亚洲欧美精品自产自拍| 日日夜夜操网爽| 男女午夜视频在线观看| 嫁个100分男人电影在线观看| 亚洲精品久久午夜乱码| 国产男女超爽视频在线观看| 国产人伦9x9x在线观看| 亚洲 国产 在线| 天堂中文最新版在线下载| 精品一区在线观看国产| 韩国精品一区二区三区| 人成视频在线观看免费观看| 亚洲av片天天在线观看| 精品人妻在线不人妻| 久久九九热精品免费| 天天躁日日躁夜夜躁夜夜| 日本欧美视频一区| 午夜福利,免费看| 九色亚洲精品在线播放| 午夜视频精品福利| 成年动漫av网址| 午夜福利免费观看在线| 免费高清在线观看视频在线观看| 国产男女内射视频| av电影中文网址| 伊人久久大香线蕉亚洲五| 人人妻人人爽人人添夜夜欢视频| 精品福利永久在线观看| 老司机午夜福利在线观看视频 | 精品一区二区三区四区五区乱码| 免费在线观看日本一区| 国产一级毛片在线| 国产男女超爽视频在线观看| 欧美精品av麻豆av| 999久久久国产精品视频| 在线看a的网站| 久久99热这里只频精品6学生| 午夜福利视频在线观看免费| 中文字幕人妻熟女乱码| 一级a爱视频在线免费观看| 捣出白浆h1v1| 国产亚洲午夜精品一区二区久久| 国产日韩欧美亚洲二区| 两性夫妻黄色片| 国产成+人综合+亚洲专区| 少妇被粗大的猛进出69影院| 国产成人欧美在线观看 | 亚洲国产精品成人久久小说| 欧美黑人欧美精品刺激| 一本—道久久a久久精品蜜桃钙片| 黑人巨大精品欧美一区二区mp4| 久久久国产精品麻豆| 日本猛色少妇xxxxx猛交久久| 午夜精品国产一区二区电影| 久久精品亚洲av国产电影网| 一本大道久久a久久精品| a级片在线免费高清观看视频| 菩萨蛮人人尽说江南好唐韦庄| 国产精品1区2区在线观看. | 天天影视国产精品| 久久性视频一级片| 色婷婷久久久亚洲欧美| 精品久久久久久久毛片微露脸 | 99久久人妻综合| 狠狠婷婷综合久久久久久88av| 精品亚洲成a人片在线观看| 久久精品aⅴ一区二区三区四区| 午夜福利在线免费观看网站| 99久久国产精品久久久| 国产1区2区3区精品| 久久久久久人人人人人| 国产精品.久久久| 亚洲国产av新网站| 精品国产国语对白av| 久久久久久免费高清国产稀缺| 久久毛片免费看一区二区三区| 性高湖久久久久久久久免费观看| 99热国产这里只有精品6| 正在播放国产对白刺激| 在线 av 中文字幕| 波多野结衣av一区二区av| 亚洲av电影在线进入| 麻豆国产av国片精品| 成人三级做爰电影| 下体分泌物呈黄色| 人妻 亚洲 视频| 麻豆av在线久日| av一本久久久久| 亚洲中文字幕日韩| 欧美黄色片欧美黄色片| 视频在线观看一区二区三区| 美女福利国产在线| 一二三四社区在线视频社区8| av网站在线播放免费| 成年av动漫网址| 亚洲精品日韩在线中文字幕| 国产精品久久久久久精品电影小说| 女人久久www免费人成看片| 男女高潮啪啪啪动态图| 精品国产一区二区三区四区第35| 国产97色在线日韩免费| 老汉色∧v一级毛片| 母亲3免费完整高清在线观看| 国产成人精品久久二区二区免费| 19禁男女啪啪无遮挡网站| av一本久久久久| 波多野结衣av一区二区av| 中文字幕人妻熟女乱码| 亚洲精品国产区一区二| 最近最新免费中文字幕在线| 五月天丁香电影| 日韩熟女老妇一区二区性免费视频| 男女高潮啪啪啪动态图| 纯流量卡能插随身wifi吗| 少妇裸体淫交视频免费看高清 | 成人三级做爰电影| 国产精品久久久久久精品电影小说| 欧美精品高潮呻吟av久久| 久久99一区二区三区| 多毛熟女@视频| 黄色 视频免费看| 国产成人精品在线电影| 欧美日韩亚洲国产一区二区在线观看 | 欧美日韩亚洲高清精品| 欧美黄色片欧美黄色片| 一级片免费观看大全| 亚洲国产精品一区三区| 91字幕亚洲| 免费高清在线观看视频在线观看| 亚洲精品一卡2卡三卡4卡5卡 | 日日爽夜夜爽网站| 99精品欧美一区二区三区四区| 亚洲avbb在线观看| 国产又爽黄色视频| 精品久久蜜臀av无| 国产欧美日韩精品亚洲av| 狂野欧美激情性xxxx| 亚洲欧美精品综合一区二区三区| 深夜精品福利| 久久亚洲精品不卡| 天天躁夜夜躁狠狠躁躁| 国产伦人伦偷精品视频| 色综合欧美亚洲国产小说| 一级黄色大片毛片| 午夜福利在线免费观看网站| 久久久久精品国产欧美久久久 | 国产又爽黄色视频| 伊人亚洲综合成人网| 国产精品一区二区在线不卡| 不卡一级毛片| 性高湖久久久久久久久免费观看| 国产成人欧美在线观看 | 老司机福利观看| 另类精品久久| 天天添夜夜摸| 国产片内射在线| 成人av一区二区三区在线看 | 脱女人内裤的视频| 91成人精品电影| 亚洲国产欧美一区二区综合| 国产高清国产精品国产三级| 亚洲专区中文字幕在线| 交换朋友夫妻互换小说| 亚洲精品av麻豆狂野| 人人妻,人人澡人人爽秒播| 成人亚洲精品一区在线观看| 美女扒开内裤让男人捅视频| 丰满饥渴人妻一区二区三| 欧美日韩亚洲综合一区二区三区_| av网站在线播放免费| 少妇粗大呻吟视频| 国产精品av久久久久免费| 国产精品亚洲av一区麻豆| 水蜜桃什么品种好| 黄色视频不卡| 少妇的丰满在线观看| 五月天丁香电影| a级片在线免费高清观看视频| 精品免费久久久久久久清纯 | 免费观看a级毛片全部| 欧美精品啪啪一区二区三区 | 国产成人免费无遮挡视频| 一二三四社区在线视频社区8| 成人18禁高潮啪啪吃奶动态图| 国产精品久久久人人做人人爽| 一级片'在线观看视频| 国产在线观看jvid| 大片电影免费在线观看免费| 久久亚洲精品不卡| 伦理电影免费视频| 他把我摸到了高潮在线观看 | 国产视频一区二区在线看| 亚洲av欧美aⅴ国产| 久久久久网色| 久久免费观看电影| 一二三四在线观看免费中文在| 伊人亚洲综合成人网| 一边摸一边抽搐一进一出视频| 男女之事视频高清在线观看| www日本在线高清视频| 男女之事视频高清在线观看| 色老头精品视频在线观看| 久久精品国产亚洲av高清一级| av电影中文网址| 亚洲九九香蕉| 天天躁日日躁夜夜躁夜夜| 黑丝袜美女国产一区| 国产成人精品久久二区二区免费| 99精国产麻豆久久婷婷| 免费在线观看视频国产中文字幕亚洲 | 青春草亚洲视频在线观看| 少妇被粗大的猛进出69影院| 精品第一国产精品| 国产欧美日韩精品亚洲av| 久久99热这里只频精品6学生| www.999成人在线观看| 色视频在线一区二区三区| 精品少妇黑人巨大在线播放| 啦啦啦在线免费观看视频4| 亚洲精品粉嫩美女一区| 精品国产国语对白av| 大陆偷拍与自拍| 中文字幕av电影在线播放| 国产av一区二区精品久久| 高清黄色对白视频在线免费看| 国产av精品麻豆| 亚洲精品国产色婷婷电影| 国产精品麻豆人妻色哟哟久久| 99久久精品国产亚洲精品| 日韩 亚洲 欧美在线| 人人妻人人澡人人爽人人夜夜| 中文字幕最新亚洲高清| 日韩视频在线欧美| 日本vs欧美在线观看视频| 少妇 在线观看| 狠狠婷婷综合久久久久久88av| 国产免费一区二区三区四区乱码| 岛国毛片在线播放| 人人澡人人妻人| 两个人免费观看高清视频| 日本撒尿小便嘘嘘汇集6| 婷婷色av中文字幕| 亚洲黑人精品在线| 丝袜美足系列| 一二三四社区在线视频社区8| 国产高清videossex| 亚洲第一av免费看| 黄色片一级片一级黄色片| 亚洲色图综合在线观看| 啦啦啦免费观看视频1| 久久精品久久久久久噜噜老黄| 亚洲五月色婷婷综合| 美女大奶头黄色视频| 老司机深夜福利视频在线观看 | 久久久久精品国产欧美久久久 | 又大又爽又粗| 一边摸一边抽搐一进一出视频| 国产有黄有色有爽视频| 免费一级毛片在线播放高清视频 | 日韩免费高清中文字幕av| 亚洲国产精品999| 欧美另类亚洲清纯唯美| 亚洲精品一区蜜桃| 免费观看a级毛片全部| 九色亚洲精品在线播放| 精品国产一区二区三区久久久樱花| 国产精品一区二区免费欧美 | 欧美大码av| 在线观看免费视频网站a站| 女人精品久久久久毛片| 乱人伦中国视频| 欧美日韩一级在线毛片| 久久久精品区二区三区| 狠狠婷婷综合久久久久久88av| 一本久久精品| 国产精品自产拍在线观看55亚洲 | 久9热在线精品视频| 亚洲欧美激情在线| 满18在线观看网站| 我要看黄色一级片免费的| 黄色视频,在线免费观看| 亚洲av成人一区二区三| 19禁男女啪啪无遮挡网站| 黄色怎么调成土黄色| 人妻 亚洲 视频| 色婷婷av一区二区三区视频| 国产高清国产精品国产三级| 国产一级毛片在线| 性高湖久久久久久久久免费观看| 亚洲欧洲日产国产| 国产亚洲精品第一综合不卡| 午夜精品久久久久久毛片777| 国产精品一区二区精品视频观看| 国产av精品麻豆| 大片免费播放器 马上看| 国产在线视频一区二区| 久久青草综合色| 亚洲精品美女久久久久99蜜臀| 亚洲国产中文字幕在线视频| 新久久久久国产一级毛片| 男人操女人黄网站| 久久久久久免费高清国产稀缺| 亚洲av日韩在线播放| 91九色精品人成在线观看| 国产不卡av网站在线观看| 少妇人妻久久综合中文| 男人爽女人下面视频在线观看| 午夜福利,免费看| 免费av中文字幕在线| 精品亚洲成国产av| 日本一区二区免费在线视频| 老司机在亚洲福利影院| 人人妻人人澡人人爽人人夜夜| 国产精品免费大片| 啦啦啦免费观看视频1| 操美女的视频在线观看| 亚洲精品国产av蜜桃| 汤姆久久久久久久影院中文字幕| a 毛片基地| 母亲3免费完整高清在线观看| 热99国产精品久久久久久7| 国产成人精品在线电影| 中文字幕av电影在线播放| 久久中文看片网| 男女国产视频网站| 亚洲成人免费av在线播放| 人人妻人人澡人人看| 国产一级毛片在线| 大香蕉久久成人网| 久久久久网色| 精品免费久久久久久久清纯 | 1024视频免费在线观看| 欧美精品啪啪一区二区三区 | 91麻豆av在线| 中文字幕制服av| 天天操日日干夜夜撸| 国产在视频线精品| 国产一卡二卡三卡精品| 国产区一区二久久| 国产成人欧美| 国产精品久久久久久人妻精品电影 | 日本wwww免费看| 欧美成人午夜精品| 肉色欧美久久久久久久蜜桃| 精品久久久精品久久久| 欧美日本中文国产一区发布| 91大片在线观看| 欧美变态另类bdsm刘玥| 91麻豆av在线| 女人高潮潮喷娇喘18禁视频| 久9热在线精品视频| 韩国高清视频一区二区三区| 电影成人av| 欧美日韩国产mv在线观看视频| 成人亚洲精品一区在线观看| 黄色毛片三级朝国网站| www.av在线官网国产| 50天的宝宝边吃奶边哭怎么回事| 宅男免费午夜| 免费久久久久久久精品成人欧美视频| tube8黄色片| 夜夜夜夜夜久久久久| 免费黄频网站在线观看国产| 一区二区av电影网| 亚洲国产欧美日韩在线播放| 黄片大片在线免费观看| 97精品久久久久久久久久精品| 黄片大片在线免费观看| av免费在线观看网站| 亚洲黑人精品在线| 免费日韩欧美在线观看| 性色av乱码一区二区三区2| 一区二区三区乱码不卡18| 叶爱在线成人免费视频播放| 亚洲av片天天在线观看| 日韩熟女老妇一区二区性免费视频| 亚洲精品国产色婷婷电影| 久久精品国产亚洲av香蕉五月 | 国产伦人伦偷精品视频| 亚洲熟女毛片儿| 免费不卡黄色视频| 亚洲精品一二三| av电影中文网址| 国产黄频视频在线观看| 日本91视频免费播放| 欧美久久黑人一区二区| 久久人妻熟女aⅴ| 欧美亚洲日本最大视频资源| 人妻久久中文字幕网| 国产欧美亚洲国产| 老司机在亚洲福利影院| 国产成人一区二区三区免费视频网站| 考比视频在线观看| 在线观看舔阴道视频| 这个男人来自地球电影免费观看| 精品亚洲成a人片在线观看| 亚洲欧洲精品一区二区精品久久久| 老司机影院毛片| 亚洲熟女毛片儿| 亚洲精品乱久久久久久| 日本黄色日本黄色录像| 老熟妇乱子伦视频在线观看 | 欧美日韩黄片免| 亚洲国产欧美在线一区| 亚洲精品国产精品久久久不卡| 老司机在亚洲福利影院| 日本wwww免费看| 宅男免费午夜| av不卡在线播放| av福利片在线| 日韩欧美一区二区三区在线观看 | 中文字幕色久视频| 久久久水蜜桃国产精品网| 欧美日韩av久久| 女警被强在线播放| tube8黄色片| 夫妻午夜视频| 一本一本久久a久久精品综合妖精| 久久人人爽av亚洲精品天堂| 久久天堂一区二区三区四区| 制服人妻中文乱码| 欧美变态另类bdsm刘玥| 久久久精品94久久精品| 日韩欧美一区视频在线观看| 国产精品欧美亚洲77777| 成年动漫av网址| 午夜成年电影在线免费观看| h视频一区二区三区| 高清av免费在线| 免费少妇av软件| 中文字幕色久视频| 美女午夜性视频免费| 免费av中文字幕在线| 考比视频在线观看| 丰满人妻熟妇乱又伦精品不卡| 欧美黑人欧美精品刺激| 中文字幕最新亚洲高清| 亚洲精品美女久久av网站| 99九九在线精品视频| 男女下面插进去视频免费观看| 建设人人有责人人尽责人人享有的| tocl精华| 精品人妻一区二区三区麻豆| 成人影院久久| 欧美黑人欧美精品刺激| 亚洲 欧美一区二区三区| 在线天堂中文资源库| 黄色毛片三级朝国网站| www.av在线官网国产| 大片免费播放器 马上看| 国产精品av久久久久免费| 99国产精品一区二区三区| 久久99一区二区三区| 咕卡用的链子| 亚洲欧美一区二区三区黑人| 黄色视频在线播放观看不卡| 男人操女人黄网站| 欧美+亚洲+日韩+国产| 十分钟在线观看高清视频www| 国产成人系列免费观看| www.精华液| 深夜精品福利| 超碰97精品在线观看|