• <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
    91精品一卡2卡3卡4卡| 国产av一区二区精品久久 | 国产片特级美女逼逼视频| 久久人人爽av亚洲精品天堂 | 狂野欧美激情性bbbbbb| 久久人人爽av亚洲精品天堂 | 免费看日本二区| 亚洲色图综合在线观看| 最近最新中文字幕免费大全7| 亚洲av二区三区四区| 小蜜桃在线观看免费完整版高清| 久热久热在线精品观看| 99久国产av精品国产电影| 国国产精品蜜臀av免费| 九色成人免费人妻av| 大香蕉久久网| 直男gayav资源| 国产成人aa在线观看| 卡戴珊不雅视频在线播放| 久久精品国产自在天天线| 黄色日韩在线| 97超碰精品成人国产| 直男gayav资源| 美女高潮的动态| 能在线免费看毛片的网站| 中文字幕精品免费在线观看视频 | 熟女av电影| 国产亚洲一区二区精品| 国精品久久久久久国模美| 草草在线视频免费看| 亚洲人成网站在线观看播放| 欧美bdsm另类| 80岁老熟妇乱子伦牲交| 国产黄色视频一区二区在线观看| 免费大片18禁| 久久鲁丝午夜福利片| 亚洲精品国产av蜜桃| 黄色怎么调成土黄色| 亚洲欧洲国产日韩| 美女脱内裤让男人舔精品视频| 亚洲精品一区蜜桃| 亚洲欧美清纯卡通| 亚洲va在线va天堂va国产| av在线播放精品| 欧美日韩视频高清一区二区三区二| 黑人高潮一二区| 99热网站在线观看| 一区二区三区精品91| 久久久久久九九精品二区国产| 美女脱内裤让男人舔精品视频| 岛国毛片在线播放| 日日摸夜夜添夜夜添av毛片| a级一级毛片免费在线观看| tube8黄色片| 直男gayav资源| 干丝袜人妻中文字幕| 国产中年淑女户外野战色| 亚洲综合精品二区| 亚洲欧美日韩卡通动漫| 日韩成人av中文字幕在线观看| 日本猛色少妇xxxxx猛交久久| 国产欧美另类精品又又久久亚洲欧美| 美女国产视频在线观看| 少妇人妻 视频| 日韩不卡一区二区三区视频在线| 免费高清在线观看视频在线观看| av在线蜜桃| 精品人妻偷拍中文字幕| 纵有疾风起免费观看全集完整版| 男女啪啪激烈高潮av片| 午夜福利高清视频| 免费人妻精品一区二区三区视频| 久久久久精品久久久久真实原创| 91aial.com中文字幕在线观看| 久久久久精品久久久久真实原创| 少妇人妻久久综合中文| 国产精品av视频在线免费观看| 最黄视频免费看| 日韩不卡一区二区三区视频在线| 少妇丰满av| 日韩精品有码人妻一区| 亚洲精品日韩av片在线观看| 国产精品久久久久久久电影| 九草在线视频观看| 中国美白少妇内射xxxbb| 亚洲三级黄色毛片| 免费黄频网站在线观看国产| 嘟嘟电影网在线观看| 18禁裸乳无遮挡免费网站照片| 18禁在线播放成人免费| 91久久精品国产一区二区三区| 九色成人免费人妻av| 日韩精品有码人妻一区| 一区在线观看完整版| 只有这里有精品99| 久久久亚洲精品成人影院| 免费观看av网站的网址| 大陆偷拍与自拍| 99国产精品免费福利视频| 久久这里有精品视频免费| 直男gayav资源| 80岁老熟妇乱子伦牲交| 亚洲三级黄色毛片| 久久久久久久国产电影| 亚洲成人一二三区av| 国产亚洲一区二区精品| a级毛色黄片| 国产精品国产三级国产专区5o| 97精品久久久久久久久久精品| 国产精品久久久久久精品电影小说 | 草草在线视频免费看| 少妇人妻一区二区三区视频| 亚洲av国产av综合av卡| 我的老师免费观看完整版| 三级国产精品欧美在线观看| 精品国产一区二区三区久久久樱花 | a级一级毛片免费在线观看| h视频一区二区三区| 久久久久久久亚洲中文字幕| 老女人水多毛片| 国产深夜福利视频在线观看| 午夜精品国产一区二区电影| 纯流量卡能插随身wifi吗| 亚洲成人一二三区av| 久热久热在线精品观看| 91精品国产九色| 女的被弄到高潮叫床怎么办| 香蕉精品网在线| 午夜日本视频在线| 久久女婷五月综合色啪小说| 国产精品一区二区在线观看99| 三级经典国产精品| av天堂中文字幕网| 国产av码专区亚洲av| 欧美日韩综合久久久久久| h日本视频在线播放| 女人久久www免费人成看片| 亚洲成色77777| 亚洲国产高清在线一区二区三| 欧美日韩一区二区视频在线观看视频在线| 久久久久久久久久人人人人人人| 少妇高潮的动态图| 欧美高清成人免费视频www| 欧美高清成人免费视频www| 亚洲无线观看免费| 精品少妇黑人巨大在线播放| 少妇高潮的动态图| 黄色配什么色好看| av福利片在线观看| 我的老师免费观看完整版| 交换朋友夫妻互换小说| 成人毛片a级毛片在线播放| av在线蜜桃| 成人毛片60女人毛片免费| 日韩人妻高清精品专区| 久久99蜜桃精品久久| 男女下面进入的视频免费午夜| 欧美日韩精品成人综合77777| 精品一区二区免费观看| av又黄又爽大尺度在线免费看| 中国三级夫妇交换| www.色视频.com| 爱豆传媒免费全集在线观看| 国产亚洲91精品色在线| 噜噜噜噜噜久久久久久91| av国产精品久久久久影院| 午夜福利视频精品| 亚洲va在线va天堂va国产| 亚洲欧美中文字幕日韩二区| 亚洲欧美日韩另类电影网站 | 熟妇人妻不卡中文字幕| 观看免费一级毛片| 亚洲人成网站高清观看| 午夜免费观看性视频| 午夜免费男女啪啪视频观看| 伦精品一区二区三区| 人人妻人人澡人人爽人人夜夜| 亚洲精品乱码久久久久久按摩| 纵有疾风起免费观看全集完整版| 肉色欧美久久久久久久蜜桃| 国产一区二区三区综合在线观看 | videos熟女内射| 伊人久久国产一区二区| 欧美成人精品欧美一级黄| 啦啦啦啦在线视频资源| 精品人妻熟女av久视频| 91精品国产国语对白视频| 成人漫画全彩无遮挡| 舔av片在线| 成人18禁高潮啪啪吃奶动态图 | 国产av国产精品国产| 视频中文字幕在线观看| 亚洲av福利一区| 国内揄拍国产精品人妻在线| 国产亚洲91精品色在线| 又粗又硬又长又爽又黄的视频| 国产一区二区在线观看日韩| av播播在线观看一区| 一区二区三区精品91| 国产精品偷伦视频观看了| 亚洲精品乱码久久久久久按摩| 久久久亚洲精品成人影院| 99精国产麻豆久久婷婷| 亚洲成色77777| 久久久久久久久久久免费av| 熟女av电影| 波野结衣二区三区在线| 国产高清三级在线| 五月开心婷婷网| 亚洲欧美一区二区三区黑人 | 日韩 亚洲 欧美在线| 久久久久久久国产电影| h视频一区二区三区| 女性生殖器流出的白浆| 成人毛片a级毛片在线播放| 日本午夜av视频| 一级a做视频免费观看| 午夜福利视频精品| 免费黄色在线免费观看| 国产无遮挡羞羞视频在线观看| 一个人免费看片子| 国产毛片在线视频| 成年美女黄网站色视频大全免费 | 欧美成人精品欧美一级黄| 精品亚洲成国产av| 内射极品少妇av片p| 国产女主播在线喷水免费视频网站| 国产一区二区在线观看日韩| 色哟哟·www| 天天躁日日操中文字幕| 日韩一区二区三区影片| 精品亚洲成国产av| 少妇的逼水好多| 亚洲av二区三区四区| 欧美亚洲 丝袜 人妻 在线| 少妇人妻久久综合中文| 天天躁日日操中文字幕| 久久精品国产a三级三级三级| 成人一区二区视频在线观看| 欧美最新免费一区二区三区| 日韩精品有码人妻一区| 激情 狠狠 欧美| 99久久综合免费| 狂野欧美白嫩少妇大欣赏| 免费黄网站久久成人精品| 国产高清三级在线| 麻豆成人av视频| 三级国产精品欧美在线观看| 3wmmmm亚洲av在线观看| 大片电影免费在线观看免费| 久久97久久精品| 日韩 亚洲 欧美在线| 亚洲美女搞黄在线观看| 欧美高清成人免费视频www| 熟女电影av网| 久久国产乱子免费精品| 精品午夜福利在线看| 久久久久久久精品精品| 国产精品熟女久久久久浪| 在线亚洲精品国产二区图片欧美 | 七月丁香在线播放| 欧美一区二区亚洲| 日韩亚洲欧美综合| 精品99又大又爽又粗少妇毛片| 男女无遮挡免费网站观看| 少妇高潮的动态图| 日日啪夜夜撸| 七月丁香在线播放| 超碰av人人做人人爽久久| 国产精品久久久久成人av| 欧美日韩亚洲高清精品| 免费人成在线观看视频色| 欧美成人午夜免费资源| 亚洲av成人精品一区久久| 内射极品少妇av片p| 欧美精品人与动牲交sv欧美| 亚洲国产av新网站| 青青草视频在线视频观看| 久久久久网色| 新久久久久国产一级毛片| 亚洲欧美日韩东京热| 日韩一本色道免费dvd| 黑人高潮一二区| 99久久综合免费| 男人狂女人下面高潮的视频| 80岁老熟妇乱子伦牲交| 高清午夜精品一区二区三区| 最黄视频免费看| 偷拍熟女少妇极品色| 国产 精品1| 99热这里只有是精品50| 亚洲色图av天堂| 免费人妻精品一区二区三区视频| 国产亚洲91精品色在线| 国模一区二区三区四区视频| 如何舔出高潮| 色5月婷婷丁香| av网站免费在线观看视频| 亚洲一级一片aⅴ在线观看| 啦啦啦视频在线资源免费观看| 日本黄色片子视频| 赤兔流量卡办理| 哪个播放器可以免费观看大片| 一级二级三级毛片免费看| 最近手机中文字幕大全| 内地一区二区视频在线| 网址你懂的国产日韩在线| 午夜福利在线观看免费完整高清在| 男女免费视频国产| 精品久久久久久久久亚洲| 亚州av有码| 久久 成人 亚洲| 一级黄片播放器| 又黄又爽又刺激的免费视频.| 久久精品国产亚洲av涩爱| 高清午夜精品一区二区三区| 蜜臀久久99精品久久宅男| 国产欧美另类精品又又久久亚洲欧美| 亚洲精品国产成人久久av| 午夜福利影视在线免费观看| 国产69精品久久久久777片| 直男gayav资源| 麻豆成人av视频| 国产免费一级a男人的天堂| 少妇人妻 视频| 舔av片在线| 国产精品嫩草影院av在线观看| 乱码一卡2卡4卡精品| 精品国产露脸久久av麻豆| 日韩三级伦理在线观看| 国产在视频线精品| 国产 一区 欧美 日韩| 夜夜看夜夜爽夜夜摸| 少妇人妻久久综合中文| 丰满少妇做爰视频| 丝袜喷水一区| 你懂的网址亚洲精品在线观看| 色吧在线观看| 女性被躁到高潮视频| 亚洲欧美中文字幕日韩二区| 大又大粗又爽又黄少妇毛片口| 亚洲精品亚洲一区二区| 亚洲婷婷狠狠爱综合网| 日日啪夜夜爽| 成人亚洲欧美一区二区av| 国产成人免费观看mmmm| 日韩免费高清中文字幕av| 丝袜喷水一区| 日本-黄色视频高清免费观看| 亚洲国产欧美在线一区| 亚洲国产精品成人久久小说| a级一级毛片免费在线观看| 亚洲精品日韩在线中文字幕| 成人18禁高潮啪啪吃奶动态图 | 久久久久久久亚洲中文字幕| 亚洲欧美精品自产自拍| 在线观看三级黄色| 国产 精品1| 亚洲精品一二三| 国产真实伦视频高清在线观看| 久久6这里有精品| 99久久人妻综合| 欧美区成人在线视频| 各种免费的搞黄视频| 国产精品嫩草影院av在线观看| 日本欧美视频一区| 久久人妻熟女aⅴ| 777米奇影视久久| 男人和女人高潮做爰伦理| 国产午夜精品一二区理论片| 哪个播放器可以免费观看大片| 亚洲va在线va天堂va国产| 成人午夜精彩视频在线观看| 国产黄色免费在线视频| 啦啦啦视频在线资源免费观看| 在线观看免费视频网站a站| 亚洲熟女精品中文字幕| 五月天丁香电影| 日韩国内少妇激情av| 美女国产视频在线观看| 国产精品国产三级国产专区5o| 老司机影院毛片| 九色成人免费人妻av| 色视频www国产| 欧美日韩综合久久久久久| 国产精品欧美亚洲77777| 亚洲国产成人一精品久久久| 最近中文字幕2019免费版| 自拍欧美九色日韩亚洲蝌蚪91 | 91久久精品国产一区二区成人| 女的被弄到高潮叫床怎么办| 国产精品99久久久久久久久| 亚洲av欧美aⅴ国产| 伦理电影免费视频| 婷婷色av中文字幕| 成人亚洲精品一区在线观看 | 啦啦啦中文免费视频观看日本| 久久久久久久国产电影| 国产熟女欧美一区二区| 高清av免费在线| 国产视频内射| 亚洲中文av在线| 爱豆传媒免费全集在线观看| 亚洲欧美中文字幕日韩二区| 国产日韩欧美亚洲二区| 中文天堂在线官网| 免费大片18禁| 国产91av在线免费观看| 日本免费在线观看一区| 夫妻午夜视频| 免费看av在线观看网站| 亚洲av二区三区四区| 久久久久视频综合| 精品一区二区三区视频在线| 亚洲欧美成人精品一区二区| 少妇被粗大猛烈的视频| 99久久精品国产国产毛片| 赤兔流量卡办理| 日本午夜av视频| 日本一二三区视频观看| 久久久久久久久久久丰满| 我的老师免费观看完整版| 欧美性感艳星| 高清毛片免费看| 九九在线视频观看精品| 国产男女内射视频| 久热久热在线精品观看| 高清午夜精品一区二区三区| 噜噜噜噜噜久久久久久91| 国产爽快片一区二区三区| 大话2 男鬼变身卡| 国产成人a∨麻豆精品| 少妇的逼水好多| 亚洲欧美日韩无卡精品| 天堂8中文在线网| 纯流量卡能插随身wifi吗| 日本-黄色视频高清免费观看| av一本久久久久| 美女视频免费永久观看网站| 欧美最新免费一区二区三区| 九色成人免费人妻av| 蜜桃亚洲精品一区二区三区| 国产精品久久久久久久电影| 伦理电影免费视频| 国产深夜福利视频在线观看| 欧美精品人与动牲交sv欧美| 91久久精品电影网| 久久久成人免费电影| 卡戴珊不雅视频在线播放| 亚洲人成网站高清观看| 成人亚洲欧美一区二区av| 国产一区亚洲一区在线观看| 啦啦啦啦在线视频资源| 亚洲成人av在线免费| 久久ye,这里只有精品| 亚洲精品一区蜜桃| 亚洲av国产av综合av卡| 国产免费一区二区三区四区乱码| 亚洲av中文字字幕乱码综合| 国产精品av视频在线免费观看| 亚洲欧美一区二区三区国产| 女性生殖器流出的白浆| 18禁裸乳无遮挡动漫免费视频| 少妇 在线观看| 国产白丝娇喘喷水9色精品| 自拍欧美九色日韩亚洲蝌蚪91 | 夜夜骑夜夜射夜夜干| 亚洲无线观看免费| 爱豆传媒免费全集在线观看| 干丝袜人妻中文字幕| 欧美xxⅹ黑人| 午夜免费鲁丝| 一级毛片aaaaaa免费看小| 在线观看一区二区三区| 狂野欧美激情性bbbbbb| 黑丝袜美女国产一区| 欧美一区二区亚洲| 在线观看一区二区三区激情| 哪个播放器可以免费观看大片| 乱系列少妇在线播放| 国产高清不卡午夜福利| 18禁在线播放成人免费| 精品一区二区三卡| 国产精品av视频在线免费观看| 韩国av在线不卡| 制服丝袜香蕉在线| 国产av国产精品国产| 国产精品人妻久久久久久| 成人国产av品久久久| 欧美日韩在线观看h| 男人爽女人下面视频在线观看| 女人久久www免费人成看片| 18禁在线播放成人免费| 欧美97在线视频| 国产精品av视频在线免费观看| 韩国av在线不卡| 国产av码专区亚洲av| 久久97久久精品| 18+在线观看网站| 黑人高潮一二区| 日日啪夜夜爽| 免费观看a级毛片全部| 国产美女午夜福利| 国产av码专区亚洲av| 亚洲色图综合在线观看| 国产成人午夜福利电影在线观看| 一本—道久久a久久精品蜜桃钙片| 亚洲精品中文字幕在线视频 | a级一级毛片免费在线观看| 亚洲伊人久久精品综合| 国产精品久久久久久精品古装| 精品人妻偷拍中文字幕| 中文字幕精品免费在线观看视频 | 大陆偷拍与自拍| 天天躁日日操中文字幕| 五月开心婷婷网| 精品久久国产蜜桃| 内地一区二区视频在线| 国产黄色免费在线视频| 亚洲va在线va天堂va国产| 欧美+日韩+精品| 大片电影免费在线观看免费| 狠狠精品人妻久久久久久综合| 91午夜精品亚洲一区二区三区| av国产久精品久网站免费入址| 18+在线观看网站| 一本—道久久a久久精品蜜桃钙片| 美女视频免费永久观看网站| 欧美精品人与动牲交sv欧美| 人人妻人人添人人爽欧美一区卜 | 岛国毛片在线播放| 国产精品爽爽va在线观看网站| 国产 一区 欧美 日韩| 3wmmmm亚洲av在线观看| 午夜激情福利司机影院| 亚洲国产高清在线一区二区三| 久久精品久久精品一区二区三区| 亚洲av中文字字幕乱码综合| 高清日韩中文字幕在线| 欧美 日韩 精品 国产| 精品人妻偷拍中文字幕| 夜夜骑夜夜射夜夜干| 少妇高潮的动态图| 青春草亚洲视频在线观看| 成人综合一区亚洲| 色婷婷av一区二区三区视频| 久久精品夜色国产| 亚洲欧美精品自产自拍| 亚洲婷婷狠狠爱综合网| 不卡视频在线观看欧美| 97在线视频观看| av在线老鸭窝| 国产伦精品一区二区三区视频9| 欧美丝袜亚洲另类| 在线看a的网站| 一区二区三区乱码不卡18| 国产成人freesex在线| 精品一区在线观看国产| 黄色欧美视频在线观看| av免费在线看不卡| 国产91av在线免费观看| 欧美日韩在线观看h| 亚洲精品第二区| 国产高潮美女av| 国产亚洲精品久久久com| 久久女婷五月综合色啪小说| 亚洲,一卡二卡三卡| 网址你懂的国产日韩在线| 日韩欧美精品免费久久| 大香蕉久久网| 视频中文字幕在线观看| av视频免费观看在线观看| 中国三级夫妇交换| 高清在线视频一区二区三区| 日日撸夜夜添| 国产色爽女视频免费观看| 国产国拍精品亚洲av在线观看| av专区在线播放| 国国产精品蜜臀av免费| 日韩欧美一区视频在线观看 | 色吧在线观看| 在线观看人妻少妇| 最近最新中文字幕免费大全7| av在线app专区| 中国三级夫妇交换| tube8黄色片| 国产男人的电影天堂91| 中文字幕久久专区| 在线观看三级黄色| 联通29元200g的流量卡| 一级片'在线观看视频| 久久精品熟女亚洲av麻豆精品| 国产精品免费大片| 一级毛片aaaaaa免费看小| 色哟哟·www| 赤兔流量卡办理| 青春草亚洲视频在线观看| av在线app专区| 国内揄拍国产精品人妻在线| 欧美高清性xxxxhd video| 夫妻午夜视频| 国产高清有码在线观看视频| 亚洲怡红院男人天堂| 免费看日本二区| 黄片wwwwww| 久久精品国产亚洲av天美| 美女内射精品一级片tv| 亚洲成人av在线免费| 婷婷色综合www| 高清黄色对白视频在线免费看 | 日韩成人av中文字幕在线观看| 99国产精品免费福利视频| 久久97久久精品| 久久久久久久亚洲中文字幕| 亚洲国产欧美在线一区|