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

    基于Models-3的自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)及其效果檢驗(yàn)

    2016-03-03 06:21:00陸成偉周來東鄧也宋丹林田紅康雪周子航胡翔
    中國環(huán)境管理 2016年2期
    關(guān)鍵詞:成都市空氣質(zhì)量修正

    陸成偉,周來東,鄧也,宋丹林,田紅,康雪,周子航,胡翔

    ( 1.成都市環(huán)境科學(xué)保護(hù)研究院, 成都 610072; 2.成都市氣象局, 成都 610072 )

    基于Models-3的自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)及其效果檢驗(yàn)

    陸成偉1,周來東1,鄧也1,宋丹林1,田紅1,康雪2,周子航1,胡翔1

    ( 1.成都市環(huán)境科學(xué)保護(hù)研究院, 成都 610072; 2.成都市氣象局, 成都 610072 )

    本文介紹了一個(gè)以Models-3為基礎(chǔ)的自動(dòng)化空氣質(zhì)量數(shù)值預(yù)報(bào)系統(tǒng),該系統(tǒng)通過Gambas、Yabasic和R語言等工具進(jìn)行開發(fā),集成WRF-SMOKE-CMAQ三個(gè)模式,可通過監(jiān)測(cè)數(shù)據(jù)進(jìn)行自動(dòng)修正,完成空氣質(zhì)量業(yè)務(wù)數(shù)值預(yù)報(bào),并將結(jié)果發(fā)布到Web服務(wù)器上進(jìn)行呈現(xiàn)。該系統(tǒng)對(duì)硬件的要求較低,將部署于一臺(tái)DELL Optiplex 9010工作站上,設(shè)置6km—2km雙層嵌套,進(jìn)行成都市空氣質(zhì)量數(shù)值預(yù)報(bào)。本文分析了成都市2014年1月1日至2014年12月31日的空氣質(zhì)量數(shù)值預(yù)報(bào)結(jié)果,評(píng)價(jià)系統(tǒng)對(duì)成都市NO2、SO2、PM10、PM2.5、O3、CO以及空氣質(zhì)量指數(shù)(AQI)的預(yù)報(bào)效果。結(jié)果顯示,系統(tǒng)對(duì)于成都市2014年空氣質(zhì)量變化情況趨勢(shì)的預(yù)報(bào)效果較好,302天有效預(yù)報(bào)中,24小時(shí)直接預(yù)報(bào)的空氣質(zhì)量等級(jí)準(zhǔn)確率為58.27%,AQI預(yù)報(bào)相關(guān)系數(shù)0.71,觀測(cè)值自動(dòng)修正預(yù)報(bào)對(duì)24小時(shí)空氣質(zhì)量預(yù)報(bào)具有明顯改善效果,使其等級(jí)預(yù)報(bào)準(zhǔn)確率達(dá)到64.9%,相關(guān)系數(shù)提高到0.89。

    成都市;Models-3;空氣質(zhì)量預(yù)報(bào);自動(dòng)修正

    引言

    成都市位于四川盆地西部,東北部為德陽市,西南與雅安市相連,南接眉山,東南毗鄰資陽市,西北緊靠阿壩自治州,東西橫距192km,南北縱距166km,面積12 390km2,屬于內(nèi)陸城市。成都市東西兩翼高差近5000m,由于地表海拔高度差異顯著,直接造成水、熱等氣候要素的空間分布不均,西部山區(qū)溫度低于東部平原地區(qū),成都市全年靜風(fēng)頻率高,相對(duì)濕度較高,容易出現(xiàn)逆溫現(xiàn)象,不利的氣象條件造成成都市易出現(xiàn)污染積累,并發(fā)生霧霾天氣,從而使成都市對(duì)空氣質(zhì)量數(shù)值預(yù)報(bào)的需求迫切。

    本文介紹了一個(gè)基于Models-3的業(yè)務(wù)化空氣質(zhì)量數(shù)值預(yù)報(bào)系統(tǒng),該系統(tǒng)以CMAQ為核心模式,使用Yabasic語言開發(fā)相應(yīng)的業(yè)務(wù)化支持軟件,實(shí)現(xiàn)空氣質(zhì)量的業(yè)務(wù)化數(shù)值預(yù)報(bào),并將系統(tǒng)部署于小型工作站上,在較低的硬件成本上實(shí)現(xiàn)空氣質(zhì)量的業(yè)務(wù)化預(yù)報(bào)。本文最后對(duì)該系統(tǒng)2014年一年的空氣質(zhì)量數(shù)值預(yù)報(bào)結(jié)果進(jìn)行分析,評(píng)價(jià)該系統(tǒng)的預(yù)報(bào)效果。

    1 系統(tǒng)構(gòu)成

    Models-3空氣質(zhì)量模擬系統(tǒng)為美國EPA于1998年提出的,以第三代空氣質(zhì)量模式CMAQ為核心,包括MM5/ WRF氣象模式和SMOKE排放清單處理模式,其應(yīng)用范圍涵蓋與空氣質(zhì)量的數(shù)值模擬相關(guān)的評(píng)價(jià)、分析和決策支持等方面,可用于NO2、SO2、顆粒物和臭氧等問題的模擬[1-9],在國內(nèi)外也有成功將CMAQ模式用于空氣質(zhì)量數(shù)值預(yù)報(bào)的案例[10,11]。由于空氣質(zhì)量數(shù)值預(yù)報(bào)在我國尚處起步階段,各地預(yù)報(bào)員能力水平差異較大,因此在進(jìn)行自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)設(shè)計(jì)時(shí)有必要開發(fā)一個(gè)基于圖形用戶界面的可視化工具,方便預(yù)報(bào)員對(duì)業(yè)務(wù)預(yù)報(bào)系統(tǒng)進(jìn)行簡(jiǎn)單的操作和必要的控制。該預(yù)報(bào)系統(tǒng)使用GAMBAS語言[12]開發(fā)系統(tǒng)監(jiān)控模塊,GAMBAS最初由Beno?t Minisini開發(fā),為L(zhǎng)inux下的一款高效快速圖形化應(yīng)用程序開發(fā)工具,具有較好的移植性。業(yè)務(wù)預(yù)報(bào)系統(tǒng)主要完成Models-3的整個(gè)運(yùn)行流程,包括調(diào)用WRF模式實(shí)現(xiàn)氣象場(chǎng)數(shù)值預(yù)報(bào)、調(diào)用SMOKE模式動(dòng)態(tài)生成排放清單,并將氣象數(shù)據(jù)和排放數(shù)據(jù)輸入CMAQ模式,進(jìn)行大氣污染物濃度數(shù)值預(yù)報(bào)。業(yè)務(wù)預(yù)報(bào)系統(tǒng)需要滿足高穩(wěn)定性、高效率和低資源占用等要求。Yabasic語言為德國計(jì)算機(jī)工程師Marc-Oliver Ihm開發(fā),具有體積小、跨平臺(tái)的特點(diǎn),可以非常便捷地對(duì)多種模式進(jìn)行交互,且穩(wěn)定性高,故用于開發(fā)業(yè)務(wù)預(yù)報(bào)系統(tǒng)。該預(yù)報(bào)系統(tǒng)還包括:監(jiān)測(cè)數(shù)據(jù)收集模塊,用于獲取準(zhǔn)實(shí)時(shí)空氣質(zhì)量監(jiān)測(cè)數(shù)據(jù);預(yù)報(bào)產(chǎn)品處理系統(tǒng),在該系統(tǒng)中生成每日所需的預(yù)報(bào)產(chǎn)品,結(jié)合實(shí)測(cè)數(shù)據(jù)進(jìn)行修正預(yù)報(bào),并對(duì)預(yù)報(bào)數(shù)據(jù)進(jìn)行管理和回顧評(píng)價(jià);預(yù)報(bào)產(chǎn)品發(fā)布系統(tǒng),用于提供Web訪問支持。系統(tǒng)目前部署于一臺(tái)DELL Optiplex 9010工作站上,系統(tǒng)配置一枚Intel I7-3770四核CPU,配置8GB內(nèi)存和4TB硬盤用于文件備份,1TB硬盤用于系統(tǒng)運(yùn)行。

    1.1 系統(tǒng)監(jiān)控模塊

    系統(tǒng)監(jiān)控模塊為一個(gè)駐留在系統(tǒng)內(nèi)存中的程序,用于檢測(cè)監(jiān)測(cè)數(shù)據(jù)收集模塊和業(yè)務(wù)預(yù)報(bào)模塊是否正常運(yùn)行,同時(shí)可以提供基本的操作和預(yù)報(bào)結(jié)果的查看。通過該系統(tǒng),可以對(duì)空氣質(zhì)量預(yù)報(bào)系統(tǒng)進(jìn)行全面監(jiān)控和簡(jiǎn)單操作,并可進(jìn)行自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)的各項(xiàng)參數(shù)設(shè)置,包括WRF、SMOKE和CMAQ的運(yùn)行參數(shù)、系統(tǒng)預(yù)報(bào)長(zhǎng)度和系統(tǒng)插件啟用情況等。該模塊的另一重要功能為定時(shí)檢測(cè)NCEP服務(wù)器上GFS全球預(yù)報(bào)系統(tǒng)某一預(yù)報(bào)時(shí)次的輸出數(shù)據(jù)是否可以下載,并在可下載時(shí)自動(dòng)開始?xì)庀髷?shù)據(jù)的下載。該系統(tǒng)目前采用GMT00時(shí)的0.25°分辨率預(yù)報(bào)場(chǎng)作為初始場(chǎng),進(jìn)行含當(dāng)日在內(nèi)的96小時(shí)的數(shù)值預(yù)報(bào),預(yù)報(bào)進(jìn)程開始于每日中午12:00,每隔20分鐘檢測(cè)一次GFS氣象數(shù)據(jù)是否可以下載,通常情況下,使用8M寬帶可以在兩小時(shí)內(nèi)完成數(shù)據(jù)下載,隨后系統(tǒng)監(jiān)控模塊將啟動(dòng)空氣質(zhì)量的數(shù)值預(yù)報(bào)流程。為了盡可能地保證系統(tǒng)運(yùn)行的穩(wěn)定性,進(jìn)入2015年后,課題組租用阿里云服務(wù)器實(shí)現(xiàn)遠(yuǎn)端前日GMT12時(shí)1°分辨率的預(yù)報(bào)數(shù)據(jù)自動(dòng)備份。系統(tǒng)監(jiān)控模塊運(yùn)行于圖形界面,操作人員可查看系統(tǒng)運(yùn)行狀態(tài),并對(duì)自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)進(jìn)行操作,其運(yùn)行界面如圖1所示。

    圖1 系統(tǒng)監(jiān)控模塊運(yùn)行界面

    1.2 業(yè)務(wù)預(yù)報(bào)系統(tǒng)

    業(yè)務(wù)預(yù)報(bào)系統(tǒng)是自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)的核心組件,該系統(tǒng)耦合WRF模式、SMOKE模式和CMAQ模式,處理模式參數(shù)修改、運(yùn)行和文件處理,依次調(diào)用三個(gè)模式完成天氣數(shù)值預(yù)報(bào)、排放清單動(dòng)態(tài)處理和空氣質(zhì)量數(shù)值預(yù)報(bào),并處理與模型嵌套運(yùn)行相關(guān)的初始條件、邊界條件生成,運(yùn)行流程如圖2所示。

    業(yè)務(wù)預(yù)報(bào)系統(tǒng)首先按預(yù)報(bào)長(zhǎng)度、模式理化參數(shù)配置等信息進(jìn)行修改,檢驗(yàn)氣象數(shù)據(jù)下載是否包括第一個(gè)時(shí)次的初始場(chǎng),以及預(yù)報(bào)長(zhǎng)度所需的最后一個(gè)時(shí)次的預(yù)報(bào)場(chǎng),并啟動(dòng)WRF模式進(jìn)行氣象數(shù)值預(yù)報(bào)?,F(xiàn)階段自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)設(shè)置進(jìn)行雙層嵌套,外層為88×88×6km的網(wǎng)格,包括四川盆地的大部分地區(qū),內(nèi)層為121×96×2km的高分辨率網(wǎng)格,包括成都市行政區(qū)劃,WRF模式的氣象參數(shù)設(shè)置參考現(xiàn)有研究成果[16],WRF模式運(yùn)行完成后,需要對(duì)WRF模式的輸出數(shù)據(jù)進(jìn)行歸檔處理,隨后調(diào)用CMAQ模式的MCIP模塊將WRF模式的輸出數(shù)據(jù)轉(zhuǎn)換為IOAPI格式。氣象模式轉(zhuǎn)換完成后,首先調(diào)用CMAQ模式的ICON模塊,從前一天的預(yù)報(bào)結(jié)果中計(jì)算本次數(shù)值預(yù)報(bào)的各污染物初始濃度,進(jìn)行熱啟動(dòng),使每日的初始濃度具有合理的空間分布,隨后業(yè)務(wù)預(yù)報(bào)系統(tǒng)調(diào)用SMOKE模式處理排放清單,進(jìn)行時(shí)間分配、空間分配和垂直分配,生成CMAQ模式所需要的netCDF排放清單。

    圖2 業(yè)務(wù)預(yù)報(bào)系統(tǒng)運(yùn)行流程[13]

    氣象數(shù)據(jù)和排放清單數(shù)據(jù)處理完成后,業(yè)務(wù)預(yù)報(bào)模塊將分別對(duì)外層網(wǎng)格Domain 1和內(nèi)層網(wǎng)格Domain 2進(jìn)行空氣質(zhì)量數(shù)值預(yù)報(bào)。處理外層網(wǎng)格預(yù)報(bào)流程時(shí),使用基于統(tǒng)計(jì)數(shù)據(jù)的網(wǎng)格邊界濃度,隨后使用CCTM模塊完成數(shù)值預(yù)報(bào)過程,并在線計(jì)算光化學(xué)反應(yīng)速率和干沉降速率。Domain1計(jì)算完成后,業(yè)務(wù)預(yù)報(bào)系統(tǒng)通過BCON模塊從Domain1的預(yù)報(bào)結(jié)果中獲取Domain 2的邊界輸入輸出濃度,進(jìn)而完成成都市轄區(qū)空氣質(zhì)量數(shù)值預(yù)報(bào)。由于CCTM模塊輸出的污染物濃度主要為體積濃度,且以模型組分為主,故需要通過模型物種計(jì)算器對(duì)污染物物種進(jìn)行合并和單位轉(zhuǎn)換,模型物種計(jì)算器基于CMAQ的COMBINE組件,除COMBINE組件中所包含的大氣污染物濃度外,還同時(shí)計(jì)算相對(duì)濕度、降水、干濕沉降、地面風(fēng)速、風(fēng)向、溫度等基本要素。

    1.3 監(jiān)測(cè)數(shù)據(jù)收集模塊

    自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)包含一個(gè)準(zhǔn)實(shí)時(shí)監(jiān)測(cè)數(shù)據(jù)收集模塊,使用Python語言和SGMLParser[14,15]開發(fā),該模塊定時(shí)解析成都市環(huán)境監(jiān)測(cè)中心站實(shí)時(shí)空氣質(zhì)量監(jiān)測(cè)數(shù)據(jù)發(fā)布網(wǎng)頁,提取相關(guān)數(shù)據(jù)后存儲(chǔ)至SQLite數(shù)據(jù)庫中,并將該數(shù)據(jù)用于空氣質(zhì)量數(shù)值預(yù)報(bào)效果的比對(duì)和修正。

    1.4 預(yù)報(bào)數(shù)據(jù)處理系統(tǒng)

    預(yù)報(bào)數(shù)據(jù)處理系統(tǒng)主要使用R語言[17]開發(fā),通過RNetCDF庫[18]處理模型物種計(jì)算器輸出的netCDF文件,計(jì)算AQI值,并制作必要的產(chǎn)品。預(yù)報(bào)數(shù)據(jù)處理系統(tǒng)根據(jù)中國環(huán)境監(jiān)測(cè)總站出具的AQI指數(shù)計(jì)算方法[19]計(jì)算預(yù)報(bào)所得的AQI指數(shù),并繪制預(yù)報(bào)所需的氣象條件變化圖、大氣污染物濃度變化圖等信息。

    預(yù)報(bào)數(shù)據(jù)處理模塊生成的產(chǎn)品包括成都市中心城區(qū)未來三天空氣質(zhì)量直接預(yù)報(bào)結(jié)果、四川盆地PM2.5逐時(shí)變化趨勢(shì)圖、各預(yù)報(bào)點(diǎn)位污染物濃度逐時(shí)變化序列圖、各預(yù)報(bào)點(diǎn)位氣象要素逐時(shí)變化序列圖、地面及四個(gè)氣壓層風(fēng)速、溫度、相對(duì)濕度等要素的逐時(shí)變化圖,以及垂直剖面上的溫度、風(fēng)速和相對(duì)濕度的逐時(shí)變化圖。此外,預(yù)報(bào)數(shù)據(jù)處理模塊還調(diào)用ARWPOST和Grads繪制地面風(fēng)場(chǎng)、地面溫度、剖面溫度、剖面風(fēng)速、剖面相對(duì)濕度以及925hPa、850hPa、700hPa、500hPa四個(gè)氣壓層上的風(fēng)場(chǎng)和濕度場(chǎng),調(diào)用Verdi程序繪制地面污染物濃度逐時(shí)空間分布,并生成GIF格式動(dòng)畫,部分產(chǎn)品如圖3所示。

    預(yù)報(bào)數(shù)據(jù)處理系統(tǒng)的另一重要功能則是實(shí)現(xiàn)空氣質(zhì)量跟蹤自動(dòng)修正預(yù)報(bào),修正預(yù)報(bào)基于現(xiàn)階段監(jiān)測(cè)數(shù)據(jù)。由于目前排放源清單統(tǒng)計(jì)數(shù)據(jù)質(zhì)量和制作水平的制約,使用數(shù)值預(yù)報(bào)提供未來空氣質(zhì)量的變化趨勢(shì)信息,并在此基礎(chǔ)上通過監(jiān)測(cè)數(shù)據(jù)進(jìn)行濃度修正,可以有效減少重污染天氣的漏報(bào),自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)在每日上午9:00至10:00間出具含當(dāng)日在內(nèi)的四日空氣質(zhì)量數(shù)值預(yù)報(bào),因此編寫程序獲取出具預(yù)報(bào)時(shí)刻前8小時(shí)的空氣質(zhì)量實(shí)測(cè)數(shù)據(jù)與直接預(yù)報(bào)中前8小時(shí)的預(yù)報(bào)濃度進(jìn)行計(jì)算,計(jì)算二者相關(guān)系數(shù)R,并分情況進(jìn)行空氣質(zhì)量數(shù)值預(yù)報(bào)修正,如下式:

    圖3 自修正空氣質(zhì)量預(yù)報(bào)空氣質(zhì)量預(yù)報(bào)產(chǎn)品

    式中,Cmod_out為修正預(yù)報(bào)修正后的預(yù)報(bào)濃度;Cmod為模式直接輸出的預(yù)報(bào)濃度;Cmod_8為預(yù)報(bào)出具的時(shí)間前8小時(shí)的逐時(shí)模式輸出濃度;Cobs_8為預(yù)報(bào)出具的時(shí)間前8小時(shí)的逐時(shí)實(shí)測(cè)濃度;計(jì)算Cmod_8和Cobs_8的相關(guān)系數(shù)R,當(dāng)相關(guān)系數(shù)R≥0.5時(shí),認(rèn)為兩組數(shù)據(jù)達(dá)到中等相關(guān)性,并在此基礎(chǔ)上求得Cmod_8和Cobs_8的一元線性擬合方程的斜率a和截距b,并使用a和b對(duì)模式輸出濃度Cmod進(jìn)行修正;當(dāng)相關(guān)系數(shù)R<0.5時(shí),認(rèn)為兩組數(shù)據(jù)相關(guān)性不佳,此時(shí)以兩組數(shù)據(jù)均值的比例計(jì)算得到濃度修正系數(shù)k,并以此對(duì)預(yù)報(bào)數(shù)據(jù)進(jìn)行修正。運(yùn)用兩種修正方法修正后的六種污染物濃度再按照AQI計(jì)算方法進(jìn)行計(jì)算,得到修正預(yù)報(bào)。

    預(yù)報(bào)數(shù)據(jù)處理系統(tǒng)每天評(píng)價(jià)前日空氣質(zhì)量直接預(yù)報(bào)的效果,該功能主要使用R語言實(shí)現(xiàn),部分統(tǒng)計(jì)和繪圖功能來自openair庫[20],使用RMarkdown庫[21]生成回顧評(píng)價(jià)報(bào)告。評(píng)價(jià)報(bào)告包括預(yù)報(bào)濃度數(shù)據(jù)對(duì)比、主要?dú)庀髤?shù)回顧、污染物濃度變化時(shí)間序列、統(tǒng)計(jì)指標(biāo)分析和預(yù)報(bào)對(duì)比,其中預(yù)報(bào)濃度數(shù)據(jù)對(duì)比列出各污染物逐時(shí)濃度變化數(shù)據(jù),主要?dú)庀髤?shù)回顧則主要給出模式預(yù)報(bào)的風(fēng)玫瑰圖、邊界層高度兩個(gè)指標(biāo),污染物濃度變化時(shí)間序列則以圖表形式顯示預(yù)報(bào)數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的變化情況,統(tǒng)計(jì)指標(biāo)分析則從FAC2系數(shù)、平均偏差、平均誤差、標(biāo)準(zhǔn)化平均偏差、標(biāo)準(zhǔn)化平均誤差、均方根誤差和相關(guān)系數(shù)R等方面對(duì)預(yù)報(bào)效果進(jìn)行綜合評(píng)價(jià),最后將預(yù)報(bào)的AQI報(bào)表與實(shí)測(cè)的AQI報(bào)表進(jìn)行對(duì)比。建立空氣質(zhì)量預(yù)報(bào)回顧評(píng)價(jià)體系有助于分析空氣質(zhì)量預(yù)報(bào)中存在的問題,以便改進(jìn)預(yù)報(bào)效果,因此回顧評(píng)價(jià)以未進(jìn)行修正的直接預(yù)報(bào)結(jié)果為評(píng)價(jià)對(duì)象。

    1.5 空氣質(zhì)量預(yù)報(bào)發(fā)布模塊

    空氣質(zhì)量預(yù)報(bào)發(fā)布模塊由三部分組成,分別為基于Linux Msmtp和Mutt構(gòu)架的電子郵件發(fā)布系統(tǒng)[22]、基于Apache建立的Web服務(wù)器[23]和使用百度地圖API的開發(fā)的WEBGIS發(fā)布系統(tǒng)。

    Linux操作系統(tǒng)下,Msmtp和Mutt的組合經(jīng)常被用于設(shè)備的自動(dòng)郵件預(yù)警,使用這兩個(gè)程序?qū)崿F(xiàn)郵件的發(fā)送系統(tǒng)不需要建立本地郵件服務(wù)器,僅需使用支持smtp認(rèn)證的公共郵箱即可,省去了本地郵件服務(wù)器的架構(gòu)和維護(hù)工作。Apache服務(wù)器用于支持通過瀏覽器對(duì)預(yù)報(bào)結(jié)果進(jìn)行訪問,配置本地路徑為Apache服務(wù)器的訪問路徑,并開放80端口用于架設(shè)Web服務(wù)器,并通過NAT地址映射實(shí)現(xiàn)互聯(lián)網(wǎng)訪問。百度地圖API[24]則提供了一套免費(fèi)的、詳細(xì)的在線地圖資源,使用百度地圖API開發(fā)的GIS發(fā)布系統(tǒng)無須授權(quán)費(fèi)用,可直觀展現(xiàn)空氣質(zhì)量數(shù)值預(yù)報(bào)的結(jié)果,百度地圖可以顯示行政區(qū)劃圖和矢量標(biāo)注,標(biāo)注物的顏色符合環(huán)境空氣質(zhì)量指數(shù)(AQI)技術(shù)規(guī)定。如圖4所示。

    圖4 基于Baidu地圖API的預(yù)報(bào)發(fā)布系統(tǒng)

    2 自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)效果評(píng)價(jià)

    自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)于2013年12月15日開始投入試運(yùn)行,2014年1月1日開始提供業(yè)務(wù)預(yù)報(bào),并業(yè)務(wù)化運(yùn)行至今,取2014年預(yù)報(bào)數(shù)據(jù),共出具有效預(yù)報(bào)302天,參照我國AQI評(píng)價(jià)標(biāo)準(zhǔn),對(duì)24小時(shí)直接預(yù)報(bào)的NO2日均值、SO2日均值、PM10日均值、PM2.5日均值、O3小時(shí)最大、O38小時(shí)滑動(dòng)最大和CO日均值進(jìn)行統(tǒng)計(jì)評(píng)價(jià),分別評(píng)價(jià)其FAC2系數(shù)、平均偏差MB、平均誤差MGE、歸一化平均偏差NMB、歸一化平均誤差NMGE、均方根誤差RMSE以及相關(guān)系數(shù)R,結(jié)果如表1所示。

    表1 不同污染物預(yù)報(bào)效果統(tǒng)計(jì)分析

    302個(gè)有效預(yù)報(bào)數(shù)據(jù)中,同時(shí)存在監(jiān)測(cè)數(shù)據(jù)的預(yù)報(bào)數(shù)量為297個(gè),因此參與統(tǒng)計(jì)的各污染物預(yù)報(bào)有效樣本個(gè)數(shù)297個(gè)。FAC2系數(shù)為預(yù)報(bào)濃度在實(shí)測(cè)濃度的0.5~2倍的數(shù)量比例,可見除SO2的FAC2系數(shù)偏低外,其余各預(yù)報(bào)指標(biāo)的FAC2系數(shù)均在80%左右,部分預(yù)報(bào)項(xiàng)目的FAC2系數(shù)達(dá)到90%以上,說明預(yù)報(bào)濃度和實(shí)測(cè)濃度的偏差范圍是可以接受的。平均偏差MB顯示,除SO2整體偏高為9.68外,其余指標(biāo)均有偏低現(xiàn)象,其中以O(shè)3小時(shí)最大濃度的偏低最為明顯,其次為PM10和PM2.5;歸一化平均偏差顯示,SO2偏高達(dá)50%以上,而PM2.5偏低達(dá)到21.08%,O3小時(shí)最大濃度的偏低幅度也達(dá)到了18.69%。對(duì)比平均誤差MGE和平均偏差MB,可見SO2指標(biāo)的MGE與|MB|差異不大,說明SO2的偏差屬于系統(tǒng)性偏高,而PM10、PM2.5、O3小時(shí)最大、O38小時(shí)滑動(dòng)最大等指標(biāo)均有明顯差異,可見這些指標(biāo)的預(yù)報(bào)同時(shí)存在偏高和偏低的現(xiàn)象。均方根誤差RMSE中,以PM10最為顯著,其次分別為O3小時(shí)最大和PM2.5,可見這些指標(biāo)預(yù)報(bào)濃度均值與對(duì)應(yīng)實(shí)測(cè)均值之間的差異波動(dòng)較大;相關(guān)系數(shù)R則表明趨勢(shì)預(yù)報(bào)最佳的指標(biāo)為PM2.5,其次為PM10和SO2。各指標(biāo)時(shí)間序列圖和散點(diǎn)圖分別如圖5至圖11所示。

    從日均值變化趨勢(shì)上看,NO2的預(yù)報(bào)日均值與對(duì)應(yīng)監(jiān)測(cè)值日均值趨勢(shì)較為符合,但在6~9月存在預(yù)報(bào)濃度偏低的現(xiàn)象,從散點(diǎn)圖上也可看出NO2的預(yù)報(bào)日均值存在偏低現(xiàn)象;SO2的預(yù)報(bào)結(jié)果中,前半年日均值變化趨勢(shì)與實(shí)測(cè)值較為接近,后半年則出現(xiàn)明顯偏高的現(xiàn)象,從散點(diǎn)分布上看,部分預(yù)報(bào)結(jié)果高于同期實(shí)測(cè)值的2倍;PM10的預(yù)報(bào)結(jié)果趨勢(shì)與監(jiān)測(cè)值的趨勢(shì)較為接近,從散點(diǎn)分布上可以看出,監(jiān)測(cè)值濃度較高的時(shí)候預(yù)報(bào)值存在一定的偏低;PM2.5的預(yù)報(bào)結(jié)果趨勢(shì)與實(shí)測(cè)日均值較為接近,和PM10一樣,在監(jiān)測(cè)值較高的時(shí)段存在一定的低報(bào)現(xiàn)象。O3兩個(gè)指標(biāo)均存在明顯的偏低現(xiàn)象,和實(shí)測(cè)數(shù)據(jù)相比,存在明顯的偏低現(xiàn)象,且并未呈現(xiàn)夏季較高的現(xiàn)象,這可能和預(yù)報(bào)系統(tǒng)垂直高度層設(shè)置有關(guān),擬在后期降低模式第一層的高度;CO的預(yù)報(bào)結(jié)果與實(shí)測(cè)結(jié)果的趨勢(shì)較為一致,但相對(duì)而言波動(dòng)較大,進(jìn)入秋季以后,預(yù)報(bào)濃度和實(shí)測(cè)濃度相比存在一定的偏高現(xiàn)象。

    分別統(tǒng)計(jì)24小時(shí)和96小時(shí)的直接預(yù)報(bào)和修正預(yù)報(bào),與實(shí)測(cè)AQI數(shù)據(jù)對(duì)比,繪制時(shí)間序列圖,如圖12所示。

    分別計(jì)算各預(yù)報(bào)AQI產(chǎn)品和實(shí)測(cè)AQI的平均偏差MB、平均誤差MGE、歸一化平均偏差NMB和歸一化平均誤差NMGE,以及均方根誤差RMSE和相關(guān)系數(shù)R,結(jié)果如表2所示。

    表2 不同預(yù)報(bào)時(shí)長(zhǎng)效果對(duì)比

    可見,受限于修正預(yù)報(bào)得出的預(yù)報(bào)結(jié)果為線性修正后的結(jié)果,對(duì)于高污染天氣和優(yōu)良天氣而言,96小時(shí)的修正預(yù)報(bào)上存在明顯的滯后現(xiàn)象,而24小時(shí)修正預(yù)報(bào)則與實(shí)測(cè)結(jié)果的趨勢(shì)較為吻合,相關(guān)系數(shù)R的結(jié)果顯示,修正預(yù)報(bào)技術(shù)對(duì)24小時(shí)預(yù)報(bào)趨勢(shì)的改善效果明顯,相關(guān)系數(shù)達(dá)到0.89,而24小時(shí)直接為0.71,同時(shí)AQI預(yù)報(bào)的平均偏差和平均誤差均得到了明顯改善,RMSE的減小也說明修正預(yù)報(bào)技術(shù)降低了24小時(shí)預(yù)報(bào)AQI與實(shí)際AQI的差異波動(dòng)。96小時(shí)預(yù)報(bào)為目前自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)的最長(zhǎng)預(yù)報(bào)時(shí)長(zhǎng),直接預(yù)報(bào)相關(guān)系數(shù)為0.66,而修正預(yù)報(bào)技術(shù)對(duì)96小時(shí)預(yù)報(bào)的趨勢(shì)產(chǎn)生了一定的影響,使其相關(guān)系數(shù)R下降為0.5,雖然修正預(yù)報(bào)有助于減小96小時(shí)預(yù)報(bào)的平均偏差,但同時(shí)也增大了其平均誤差和均方根誤差,導(dǎo)致96小時(shí)修正預(yù)報(bào)與實(shí)際結(jié)果間差異波動(dòng)增大。

    3 結(jié)論

    (1)基于Models-3開發(fā)的自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)可以在較低硬件投入下實(shí)現(xiàn)空氣質(zhì)量的業(yè)務(wù)化數(shù)值預(yù)報(bào),該系統(tǒng)可以直觀進(jìn)行模擬參數(shù)的配置,具備較高的可系統(tǒng)移植,插件系統(tǒng)在不影響主模塊功能的情況下對(duì)系統(tǒng)功能進(jìn)行擴(kuò)展,具有較高的靈活性。

    (2)2014年,24小時(shí)直接預(yù)報(bào)等級(jí)準(zhǔn)確天數(shù)為176天,等級(jí)和首要污染物均準(zhǔn)確的天數(shù)為63天,24小時(shí)修正預(yù)報(bào)等級(jí)準(zhǔn)確天數(shù)為196天,等級(jí)和首要污染物均準(zhǔn)確的天數(shù)為121天,而96小時(shí)直接預(yù)報(bào)等級(jí)準(zhǔn)確天數(shù)為167天,等級(jí)和首要污染物均準(zhǔn)確的天數(shù)為58天,96小時(shí)修正預(yù)報(bào)等級(jí)準(zhǔn)確天數(shù)為112天,等級(jí)和首要污染物均準(zhǔn)確的天數(shù)為47天,可見使用實(shí)測(cè)數(shù)據(jù)對(duì)24小時(shí)數(shù)值預(yù)報(bào)進(jìn)行修正后可明顯改善其預(yù)報(bào)效果,但本文所用的修正方法對(duì)96小時(shí)預(yù)報(bào)的預(yù)報(bào)效果無明顯改善效果。

    (3)目前采用的預(yù)報(bào)修正技術(shù)存在一定的局限性,下一步工作中應(yīng)探索模式預(yù)報(bào)的氣象、空氣質(zhì)量數(shù)據(jù)和相關(guān)實(shí)測(cè)值之間的關(guān)系,建立更為科學(xué)的預(yù)報(bào)修正技術(shù)。

    (4)系統(tǒng)可以在一定程度上把握成都市空氣質(zhì)量變化的節(jié)點(diǎn),但進(jìn)入冬季后,自修正空氣質(zhì)量預(yù)報(bào)系統(tǒng)對(duì)區(qū)域性重污染天氣的把握能力有限,可在下一步工作中增大預(yù)報(bào)范圍,引入?yún)^(qū)域清單,加強(qiáng)區(qū)域性污染天氣的預(yù)報(bào)能力。

    圖5 NO2日均值預(yù)報(bào)——實(shí)測(cè)時(shí)間序列圖和散點(diǎn)圖

    圖6 SO2日均值預(yù)報(bào)——實(shí)測(cè)時(shí)間序列圖和散點(diǎn)圖

    圖7 PM10日均值預(yù)報(bào)——實(shí)測(cè)時(shí)間序列圖和散點(diǎn)圖

    圖8 PM2.5日均值預(yù)報(bào)——實(shí)測(cè)時(shí)間序列圖和散點(diǎn)圖

    圖9 O3小時(shí)最大預(yù)報(bào)——實(shí)測(cè)時(shí)間序列圖和散點(diǎn)圖

    圖10 O38小時(shí)滑動(dòng)最大值預(yù)報(bào)——實(shí)測(cè)時(shí)間序列圖和散點(diǎn)圖

    圖11 CO預(yù)報(bào)——實(shí)測(cè)日均值時(shí)間序列圖和散點(diǎn)圖

    [1] 翟世賢, 安興琴, 劉俊, 等. 不同時(shí)刻污染減排對(duì)北京市PM2.5濃度的影響[J]. 中國環(huán)境科學(xué), 2014, 34(06): 1369-1379.

    [2] Shrestha K L, Konda A, Kaga A, et al. High-resolution modeling and evaluation of ozone air quality of Osaka using MM5-CMAQ system[J]. Journal of environmental sciences, 2009, 21(06): 782-789.

    [3] 劉俊, 安興琴, 朱彤, 等. 京津冀及周邊減排對(duì)北京市PM2.5濃度下降評(píng)估研究[J]. 中國環(huán)境科學(xué), 2014, 34(11): 2726-2733.

    [4] 張艷, 余綺, 伏晴艷, 等. 長(zhǎng)江三角洲區(qū)域輸送對(duì)上海市空氣質(zhì)量影響的特征分析[J]. 中國環(huán)境科學(xué), 2010, 30(07): 914-923.

    [5] 權(quán)建農(nóng), 張曉山, 張薔, 等. 我國燃煤汞沉降的數(shù)值模擬[J]. 高原氣象, 2009, 28(01): 159-164.

    圖12 24小時(shí)和96小時(shí)直接預(yù)報(bào)與修正預(yù)報(bào)效果對(duì)比

    [6] MAI K, RYOZO O al. A numerical study of summer ozone concentration over the Kanto area of Japan using the MM5/ CMAQ model[J]. Journal of environmental sciences, 2011, 23(02): 236-246.

    [7] ZHANG M G, HAN Z W, ZHU L Y. Simulation of atmospheric aerosols in East Asia using modeling system RAMS-CMAQ: model evaluation[J]. China particuology, 2007, 5(05): 321-327.

    [8] JOSE R S , PERE J L, MORANT J L, et al. The use of modern third-generation air quality models (MM5-EMIMO-CMAQ) for real-time operational air quality impact assessment of industrial plants[J]. Water, Air, & Soil Pollution: Focus, 2009, 9(1-2): 27-37.

    [9] Shi C E, Zhang B N. Tropospheric NO2columns over Northeastern North America: comparison of CMAQ model simulations with GOME satellite measurements[J]. Advances in Atmospheric Sciences, 2008, 25(01): 59-71.

    [10] 陳彬彬, 林長(zhǎng)城, 楊凱, 等. 基于CMAQ模式產(chǎn)品的福州市空氣質(zhì)量預(yù)報(bào)系統(tǒng)[J]. 中國環(huán)境科學(xué), 2012, 32(10): 1744-1752.

    [11] 許建明, 徐祥德, 劉煜, 等. CMAQ-MOS區(qū)域空氣質(zhì)量統(tǒng)計(jì)修正模型預(yù)報(bào)途徑研究[J]. 中國科學(xué)D輯: 地球科學(xué), 2005, (35)(SI): 131-144.

    [12] Bain M A. An introduction to Gambas[J]. Linux Journal, 2006, 2006(146): 1-5.

    [13] USEPA. Operational Guidance for the Community Multiscale Air Quality (CMAQ) Modeling System. [M/OL]. Chapel Hill, NC: University of North Carolina at Chapel Hill, 2010. [2010-09-11].http://www.cmascenter.org/cmaq/documentation /4.7.1/operational-Guidance-Document.pdf.

    [14] Python software foundation. SGMLlib Simple SGML parser[EB/OL].[2012-04-21]. https://docs.python.org/2/ library/sgmllib.html#module-sgmllib.

    [15] 王琳琳. 基于HTML Parser的Web信息提取技術(shù)[D]. 北京: 北京郵電大學(xué), 2007.

    [16] 姚琳, 葉芝祥, 陸成偉, 等. 成都市空氣質(zhì)量預(yù)報(bào)中WRF的本地化參數(shù)選取[J]. 成都信息工程學(xué)院學(xué)報(bào), 2012, 27(05): 485-489.

    [17] Team R C. R: A language and environment for statistical computing [M/OL]. Vienna, Austria: 2014. [2016-03-10]. https://cram.r-profect.org/doc/manuals/fullrefman.pdf.

    [18] Michna P. RNetCDF: R Interface to NetCDF Datasets. R package version 1.6. 3-1[EB/OL].[2012-08-02]. http:// CRAN.R-project.org/package=RNetCDF.

    [19] 中國環(huán)境監(jiān)測(cè)總站. 環(huán)境空氣質(zhì)量指數(shù)(AQI)技術(shù)規(guī)定[EB/OL]. [2013-11-13]. www.cnemc.cn/publish/totalWebSite/ news/news_38845.html.

    [20] Carslaw D C, Ropkins K. Openair—an R package for air quality data analysis. Environmental Modelling & Software, 2012, 27-28: 52-61.

    [21] Allaire J J, Cheng J, Xie Y H, et al. Rmarkdown: dynamic Documents for R. 2015.[EB/OL].[2014-06-03]. http:// rmarkdewn.rstudio.com/.

    [22] 李雪白. 配置mutt做郵件客服端[EB/OL]. [2011-05-30]. http://home.ustc.edu.cn/~lixuebai/GNU/MuttConfig.html.

    [23] 謝輝. 基于Linux系統(tǒng)下Apache軟件的Web服務(wù)器設(shè)置[J]. 計(jì)算機(jī)光盤軟件與應(yīng)用, 2014, (14): 78-79.

    [24] 百度地圖. Javascript API大眾版[EB/OL]. [2014-02-23]. http://bsyun.baidu.com/index.php?title=jspopular.

    A Models-3 Based Self-correcting Air Quality Forecast System and the Estimation

    Lu Chengwei1, Zhou Laidong1, Deng Ye1, Song Danlin1, Tian Hong1, Kang Xue2, Zhou Zihang1, Hu Xiang1
    ( 1.Chengdu Academy of Environmental Sciences, Chengdu 610072;2.Chengdu Meteorological Bureau, Chengdu 610072 )

    A Models-3 based self-correcting air quality forecast system was discussed in this paper, the automated system was developed with Gambas, Yabasic and R language, which integrates 3 models including WRF, SMOKE and CMAQ. The forecast system captures monitor data from network, corrects the concentrations of different pollutants, and then public the results via web server. The hardware requirements of forecast system is relatively low and it was deployed on a DELL Optiplex 9010 workstation with a 6km-2km nested domain, giving operational air quality forecast for Chengdu. A estimation of the system was performed with 2014 forecasted concentrations and AQI, the results showed that the system well reflected the air quality variations in 2014, the hit rate of 24h direct forecast on air quality grads was 58.27% with a correlation coefficient of 0.71, and the corrected 24h forecast had a hit rate of 64.9% with a correlation coefficient of 0.89, the selfcorrecting method can improve the 24h forecast of Chengdu.

    Chengdu City; Models-3; Air quality forecast; CMAQ; Auto Correction

    X831

    1674-6252(2016)02-0102-08

    A

    10.16868/j.cnki.1674-6252.2016.02.102

    成環(huán)科研ky2013第020號(hào)。

    陸成偉(1987—),助理工程師,主要研究方向?yàn)榭諝赓|(zhì)量數(shù)值模擬與預(yù)報(bào)。

    猜你喜歡
    成都市空氣質(zhì)量修正
    中共成都市第十四屆委員會(huì)常委簡(jiǎn)歷
    先鋒(2022年4期)2022-05-07 20:26:31
    Some new thoughts of definitions of terms of sedimentary facies: Based on Miall's paper(1985)
    修正這一天
    快樂語文(2021年35期)2022-01-18 06:05:30
    成都市青羊區(qū):推行“一網(wǎng)通辦”下的“最多跑一次”
    合同解釋、合同補(bǔ)充與合同修正
    法律方法(2019年4期)2019-11-16 01:07:28
    2019年1~6月成都市經(jīng)濟(jì)運(yùn)行情況
    先鋒(2019年8期)2019-09-09 06:35:59
    2018年1—12月成都市經(jīng)濟(jì)運(yùn)行情況
    先鋒(2019年2期)2019-03-27 09:31:22
    軟件修正
    “空氣質(zhì)量發(fā)布”APP上線
    車內(nèi)空氣質(zhì)量標(biāo)準(zhǔn)進(jìn)展
    汽車與安全(2016年5期)2016-12-01 05:22:14
    天堂网av新在线| 亚洲图色成人| 国产v大片淫在线免费观看| 少妇人妻一区二区三区视频| 久久久久久久久中文| 伦理电影大哥的女人| 亚洲中文字幕一区二区三区有码在线看| 99国产极品粉嫩在线观看| 免费av毛片视频| 久久欧美精品欧美久久欧美| ponron亚洲| 亚洲精品自拍成人| 身体一侧抽搐| 日本成人三级电影网站| 亚洲国产高清在线一区二区三| 少妇熟女aⅴ在线视频| 精品99又大又爽又粗少妇毛片| 亚洲一区高清亚洲精品| 亚洲欧洲日产国产| 久久久久久国产a免费观看| 色播亚洲综合网| 青青草视频在线视频观看| 一边摸一边抽搐一进一小说| 男人的好看免费观看在线视频| 可以在线观看的亚洲视频| 国产亚洲欧美98| 老女人水多毛片| 人人妻人人看人人澡| 国产伦理片在线播放av一区 | 国产片特级美女逼逼视频| 久久精品国产亚洲网站| 如何舔出高潮| av专区在线播放| 少妇熟女aⅴ在线视频| 国产真实伦视频高清在线观看| 久久久精品欧美日韩精品| 亚洲国产色片| 日韩一区二区视频免费看| 免费观看在线日韩| 爱豆传媒免费全集在线观看| 亚洲av熟女| 最后的刺客免费高清国语| 麻豆一二三区av精品| 亚洲国产精品成人久久小说 | 蜜桃久久精品国产亚洲av| 免费搜索国产男女视频| 男插女下体视频免费在线播放| 久久精品综合一区二区三区| 久久久久久久久久成人| 在线观看美女被高潮喷水网站| 99热6这里只有精品| 中文精品一卡2卡3卡4更新| 欧美成人精品欧美一级黄| 成人鲁丝片一二三区免费| 看片在线看免费视频| 人人妻人人澡人人爽人人夜夜 | 久久久久免费精品人妻一区二区| 麻豆成人av视频| 少妇熟女aⅴ在线视频| 亚洲电影在线观看av| 久久久久久久久久成人| 欧美日韩综合久久久久久| 精品久久久久久久人妻蜜臀av| 天堂网av新在线| 51国产日韩欧美| 91久久精品国产一区二区成人| АⅤ资源中文在线天堂| 美女大奶头视频| 国产女主播在线喷水免费视频网站 | 国产69精品久久久久777片| 99久久中文字幕三级久久日本| 精品久久国产蜜桃| 久久午夜福利片| 国产中年淑女户外野战色| 国产精品麻豆人妻色哟哟久久 | 六月丁香七月| 成人鲁丝片一二三区免费| 99久久精品一区二区三区| 在线a可以看的网站| 日韩av不卡免费在线播放| 婷婷精品国产亚洲av| 人人妻人人看人人澡| 欧美激情久久久久久爽电影| 亚洲第一区二区三区不卡| av国产免费在线观看| 色噜噜av男人的天堂激情| 欧美精品国产亚洲| av天堂在线播放| 欧美在线一区亚洲| 久久精品久久久久久噜噜老黄 | 免费观看a级毛片全部| 欧美日韩乱码在线| 国产精品蜜桃在线观看 | 欧美性猛交黑人性爽| 午夜激情福利司机影院| 麻豆成人av视频| 97人妻精品一区二区三区麻豆| 午夜亚洲福利在线播放| 国产v大片淫在线免费观看| 国产毛片a区久久久久| 男人舔女人下体高潮全视频| 一边摸一边抽搐一进一小说| 精品久久久噜噜| 悠悠久久av| 国产一区二区在线av高清观看| 亚洲欧美日韩高清专用| 国产精品野战在线观看| 精品一区二区三区人妻视频| 国国产精品蜜臀av免费| 日韩亚洲欧美综合| 看片在线看免费视频| 久久99精品国语久久久| 听说在线观看完整版免费高清| 久久久久性生活片| 久久久久国产网址| 国产精品久久久久久精品电影| 亚洲欧美清纯卡通| 天堂网av新在线| 久久热精品热| 舔av片在线| 非洲黑人性xxxx精品又粗又长| a级一级毛片免费在线观看| 中文字幕av成人在线电影| 亚洲精品国产成人久久av| 国产视频内射| 内射极品少妇av片p| 日韩欧美精品v在线| 久久久精品欧美日韩精品| 国产一级毛片在线| 亚洲经典国产精华液单| av又黄又爽大尺度在线免费看 | 老司机福利观看| .国产精品久久| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 看片在线看免费视频| 最好的美女福利视频网| 免费搜索国产男女视频| 精品久久久久久久久久免费视频| ponron亚洲| 欧美性感艳星| 国产91av在线免费观看| 欧美三级亚洲精品| 热99在线观看视频| 亚洲无线在线观看| 亚洲激情五月婷婷啪啪| 久久久国产成人精品二区| 久久99蜜桃精品久久| 亚洲国产色片| 国产黄色视频一区二区在线观看 | 色综合色国产| 长腿黑丝高跟| 欧美不卡视频在线免费观看| 欧美最黄视频在线播放免费| 亚洲乱码一区二区免费版| 欧美变态另类bdsm刘玥| 日本熟妇午夜| 伊人久久精品亚洲午夜| 日韩成人伦理影院| 国产精品野战在线观看| .国产精品久久| av免费观看日本| 久久久久久久午夜电影| 一本久久中文字幕| 蜜桃久久精品国产亚洲av| 日韩成人伦理影院| 欧美+亚洲+日韩+国产| 免费看a级黄色片| 国产美女午夜福利| 日韩大尺度精品在线看网址| 看黄色毛片网站| 边亲边吃奶的免费视频| av在线观看视频网站免费| 欧美一区二区国产精品久久精品| 天天躁日日操中文字幕| 日本五十路高清| 亚洲欧美日韩东京热| 久久6这里有精品| 日韩欧美精品免费久久| 人妻系列 视频| 麻豆一二三区av精品| 男的添女的下面高潮视频| 高清日韩中文字幕在线| 日韩高清综合在线| 亚洲图色成人| 岛国毛片在线播放| 哪个播放器可以免费观看大片| 舔av片在线| 一级毛片我不卡| 一级毛片aaaaaa免费看小| 亚洲精品自拍成人| 97人妻精品一区二区三区麻豆| 高清毛片免费看| 久久久久久久亚洲中文字幕| 久久久久免费精品人妻一区二区| 熟女电影av网| 亚洲自偷自拍三级| 99热全是精品| 三级男女做爰猛烈吃奶摸视频| 国产成人精品婷婷| 精品99又大又爽又粗少妇毛片| 人妻夜夜爽99麻豆av| 校园人妻丝袜中文字幕| 精品国内亚洲2022精品成人| eeuss影院久久| 麻豆国产av国片精品| 天天一区二区日本电影三级| 亚洲内射少妇av| 婷婷色av中文字幕| 寂寞人妻少妇视频99o| 成人特级av手机在线观看| 99热这里只有精品一区| 午夜福利高清视频| 六月丁香七月| 尾随美女入室| 国产高清视频在线观看网站| 91精品国产九色| 日本-黄色视频高清免费观看| 青春草亚洲视频在线观看| 精品久久久久久久久久久久久| 亚洲精品日韩在线中文字幕 | 观看免费一级毛片| 亚洲丝袜综合中文字幕| 日韩三级伦理在线观看| www.色视频.com| 久久久精品欧美日韩精品| 美女cb高潮喷水在线观看| 蜜桃久久精品国产亚洲av| 免费观看a级毛片全部| 国产精品一区二区在线观看99 | 亚洲av第一区精品v没综合| 久久精品人妻少妇| 蜜臀久久99精品久久宅男| 日韩精品有码人妻一区| 我要看日韩黄色一级片| 日韩一区二区三区影片| 亚洲精华国产精华液的使用体验 | 99国产精品一区二区蜜桃av| 极品教师在线视频| 97超碰精品成人国产| 亚洲综合色惰| a级毛色黄片| 熟女电影av网| 欧美zozozo另类| 中文字幕av在线有码专区| 亚洲av.av天堂| 最近的中文字幕免费完整| 在线a可以看的网站| 欧美一区二区亚洲| 天堂网av新在线| 久久精品久久久久久噜噜老黄 | 18+在线观看网站| 久久99热这里只有精品18| 三级毛片av免费| 婷婷亚洲欧美| 成人漫画全彩无遮挡| 老师上课跳d突然被开到最大视频| 欧美日韩在线观看h| av女优亚洲男人天堂| 久久久久久久久久久免费av| av在线天堂中文字幕| 晚上一个人看的免费电影| 91午夜精品亚洲一区二区三区| 美女xxoo啪啪120秒动态图| 亚洲在线自拍视频| 精品一区二区免费观看| 国产精品久久久久久久久免| 午夜福利高清视频| 禁无遮挡网站| 白带黄色成豆腐渣| 日韩成人av中文字幕在线观看| 免费人成在线观看视频色| 菩萨蛮人人尽说江南好唐韦庄 | 中文字幕制服av| 亚洲成人久久性| 嫩草影院入口| 国产精品99久久久久久久久| 亚洲精品国产成人久久av| 又黄又爽又刺激的免费视频.| 国产不卡一卡二| 秋霞在线观看毛片| 欧美三级亚洲精品| videossex国产| 日韩成人伦理影院| 一进一出抽搐gif免费好疼| 日日撸夜夜添| 可以在线观看毛片的网站| 亚洲成人av在线免费| 国产精品一区二区三区四区久久| 六月丁香七月| 99久国产av精品| 亚洲欧美成人精品一区二区| 亚洲在线观看片| 综合色av麻豆| 日韩av在线大香蕉| 亚洲图色成人| 成人亚洲欧美一区二区av| 亚洲欧洲国产日韩| 一区二区三区免费毛片| 九九热线精品视视频播放| 日韩三级伦理在线观看| 观看免费一级毛片| av天堂中文字幕网| 亚洲精品自拍成人| 日韩三级伦理在线观看| 两个人视频免费观看高清| 中文字幕人妻熟人妻熟丝袜美| 91久久精品电影网| a级毛片a级免费在线| 精品日产1卡2卡| 欧美日韩一区二区视频在线观看视频在线 | 美女被艹到高潮喷水动态| 亚洲成人av在线免费| 狠狠狠狠99中文字幕| 一级毛片aaaaaa免费看小| 日韩中字成人| av天堂在线播放| 国产成人一区二区在线| 变态另类成人亚洲欧美熟女| 久久精品国产亚洲av香蕉五月| 国产片特级美女逼逼视频| 91狼人影院| 久久婷婷人人爽人人干人人爱| 成人av在线播放网站| 一个人看视频在线观看www免费| 精品国产三级普通话版| 成人av在线播放网站| 岛国毛片在线播放| 国产男人的电影天堂91| 亚洲成人久久爱视频| 国产一级毛片七仙女欲春2| 一区二区三区免费毛片| 欧美最新免费一区二区三区| 美女内射精品一级片tv| 一个人看视频在线观看www免费| 乱人视频在线观看| 国产伦理片在线播放av一区 | 成人鲁丝片一二三区免费| 此物有八面人人有两片| 简卡轻食公司| 在线观看av片永久免费下载| 可以在线观看的亚洲视频| 国产女主播在线喷水免费视频网站 | 热99在线观看视频| 别揉我奶头 嗯啊视频| 美女大奶头视频| 啦啦啦啦在线视频资源| 久久久午夜欧美精品| 一进一出抽搐动态| 波多野结衣高清无吗| 最近的中文字幕免费完整| 国内少妇人妻偷人精品xxx网站| 麻豆国产av国片精品| 国产成人精品久久久久久| 久久久久九九精品影院| 日本成人三级电影网站| 亚洲欧美清纯卡通| 亚洲中文字幕一区二区三区有码在线看| 一本久久中文字幕| 永久网站在线| 在现免费观看毛片| 悠悠久久av| 久久久国产成人免费| 国产午夜福利久久久久久| 99久久精品国产国产毛片| 中文亚洲av片在线观看爽| 精品久久久久久久末码| 在线天堂最新版资源| 欧美成人一区二区免费高清观看| 欧美日韩综合久久久久久| 国产精品久久视频播放| 精品不卡国产一区二区三区| 禁无遮挡网站| 悠悠久久av| 国产精品一区www在线观看| 日韩欧美 国产精品| 久久午夜亚洲精品久久| 中文亚洲av片在线观看爽| 高清毛片免费观看视频网站| 不卡一级毛片| 国产人妻一区二区三区在| 熟女电影av网| 老司机影院成人| 亚洲欧美精品自产自拍| 亚洲欧美中文字幕日韩二区| 久久欧美精品欧美久久欧美| 女人被狂操c到高潮| 精品人妻视频免费看| 国产精品蜜桃在线观看 | 日本撒尿小便嘘嘘汇集6| 最好的美女福利视频网| 晚上一个人看的免费电影| 国产不卡一卡二| 高清毛片免费看| www.色视频.com| 最近最新中文字幕大全电影3| 国产精品免费一区二区三区在线| 联通29元200g的流量卡| 亚洲精品国产成人久久av| 国产亚洲精品av在线| av在线天堂中文字幕| 久久久久久久久久成人| 夜夜爽天天搞| 九九在线视频观看精品| 在线免费十八禁| 日本免费一区二区三区高清不卡| 99热6这里只有精品| 国产av在哪里看| 国产单亲对白刺激| 草草在线视频免费看| 中文字幕人妻熟人妻熟丝袜美| 99热只有精品国产| 成人午夜精彩视频在线观看| 综合色丁香网| 人人妻人人澡人人爽人人夜夜 | 免费不卡的大黄色大毛片视频在线观看 | 18禁在线无遮挡免费观看视频| 女同久久另类99精品国产91| 久久久国产成人免费| 91午夜精品亚洲一区二区三区| 青春草国产在线视频 | 久久午夜亚洲精品久久| 亚洲精品亚洲一区二区| 女人被狂操c到高潮| 欧美bdsm另类| 国产亚洲av片在线观看秒播厂 | 我的老师免费观看完整版| 麻豆国产av国片精品| 久久精品人妻少妇| 免费看日本二区| 久久亚洲精品不卡| 天堂网av新在线| 久久久久久久午夜电影| 天天一区二区日本电影三级| 日韩欧美精品免费久久| 国产爱豆传媒在线观看| 成年免费大片在线观看| 欧美bdsm另类| 又黄又爽又刺激的免费视频.| 久久久久久久久大av| 一级av片app| 99热这里只有是精品50| 啦啦啦观看免费观看视频高清| 亚洲成人精品中文字幕电影| 熟女人妻精品中文字幕| av天堂中文字幕网| 最近中文字幕高清免费大全6| 国产精品电影一区二区三区| 亚洲精品乱码久久久久久按摩| 少妇猛男粗大的猛烈进出视频 | 久久久a久久爽久久v久久| 久久精品国产自在天天线| 麻豆精品久久久久久蜜桃| 18禁黄网站禁片免费观看直播| 免费看美女性在线毛片视频| 青春草视频在线免费观看| 久久久久久伊人网av| 边亲边吃奶的免费视频| 三级毛片av免费| 在线观看免费视频日本深夜| 91在线精品国自产拍蜜月| 色综合色国产| 久久草成人影院| 白带黄色成豆腐渣| 成人特级av手机在线观看| 在线观看午夜福利视频| 中文在线观看免费www的网站| 久久久久久九九精品二区国产| 麻豆久久精品国产亚洲av| 日本爱情动作片www.在线观看| 91精品一卡2卡3卡4卡| av免费在线看不卡| 国产亚洲欧美98| 日本-黄色视频高清免费观看| 一区二区三区四区激情视频 | 岛国毛片在线播放| 亚洲欧美日韩东京热| 日韩制服骚丝袜av| 最近的中文字幕免费完整| 国产白丝娇喘喷水9色精品| 欧美性猛交╳xxx乱大交人| 美女被艹到高潮喷水动态| 99久久精品一区二区三区| 成人午夜精彩视频在线观看| 五月伊人婷婷丁香| 欧美色欧美亚洲另类二区| 九草在线视频观看| 女同久久另类99精品国产91| 亚洲av一区综合| 国产精品野战在线观看| 美女黄网站色视频| 级片在线观看| 精品国产三级普通话版| 午夜福利成人在线免费观看| 三级男女做爰猛烈吃奶摸视频| www日本黄色视频网| 特大巨黑吊av在线直播| 精品少妇黑人巨大在线播放 | 久久亚洲精品不卡| 51国产日韩欧美| 中文字幕av成人在线电影| 亚洲欧美日韩高清专用| 一区二区三区高清视频在线| 国产成人影院久久av| 看片在线看免费视频| 九色成人免费人妻av| 国产色婷婷99| 精品午夜福利在线看| 日本一本二区三区精品| 免费av毛片视频| 国产精品久久久久久久久免| av国产免费在线观看| 中文字幕制服av| 一区福利在线观看| 亚洲av电影不卡..在线观看| 国内揄拍国产精品人妻在线| 中文字幕免费在线视频6| 久久精品国产自在天天线| 欧美高清性xxxxhd video| 99久久精品热视频| 小蜜桃在线观看免费完整版高清| 亚洲av熟女| 成人毛片60女人毛片免费| 嘟嘟电影网在线观看| 亚洲18禁久久av| 精品少妇黑人巨大在线播放 | 亚洲人成网站高清观看| 观看美女的网站| 国产午夜福利久久久久久| 中文字幕精品亚洲无线码一区| 一区福利在线观看| 欧美变态另类bdsm刘玥| 精品久久久久久久末码| 日日啪夜夜撸| 亚洲天堂国产精品一区在线| 综合色av麻豆| 在线观看免费视频日本深夜| 日本撒尿小便嘘嘘汇集6| 亚洲av.av天堂| 日日摸夜夜添夜夜添av毛片| 毛片一级片免费看久久久久| 91午夜精品亚洲一区二区三区| 日本三级黄在线观看| 亚洲无线在线观看| 国产老妇女一区| 99热这里只有是精品在线观看| 国产精品一及| 欧美激情国产日韩精品一区| 久久午夜亚洲精品久久| eeuss影院久久| 国产色爽女视频免费观看| 国产黄a三级三级三级人| 熟妇人妻久久中文字幕3abv| 春色校园在线视频观看| 日本色播在线视频| 国产日韩欧美在线精品| 成人午夜高清在线视频| 日韩人妻高清精品专区| 国产av一区在线观看免费| 欧美性猛交╳xxx乱大交人| 夜夜看夜夜爽夜夜摸| 国产单亲对白刺激| 天堂√8在线中文| 伦理电影大哥的女人| 少妇人妻精品综合一区二区 | 精品日产1卡2卡| 麻豆国产97在线/欧美| 亚洲av二区三区四区| 99热全是精品| 99热只有精品国产| 人妻制服诱惑在线中文字幕| 国产美女午夜福利| 亚洲av.av天堂| 熟妇人妻久久中文字幕3abv| 精品一区二区三区视频在线| 亚洲欧洲国产日韩| 欧美日本视频| 久久九九热精品免费| 2021天堂中文幕一二区在线观| 最近最新中文字幕大全电影3| 精品免费久久久久久久清纯| 两个人的视频大全免费| 亚洲美女搞黄在线观看| 91久久精品国产一区二区三区| 亚洲人成网站在线播放欧美日韩| 麻豆精品久久久久久蜜桃| 黑人高潮一二区| 国产精品一区二区三区四区免费观看| 一进一出抽搐动态| 久久99蜜桃精品久久| 嘟嘟电影网在线观看| 日本黄大片高清| 国产亚洲91精品色在线| 日韩视频在线欧美| 亚洲国产欧美人成| 久久久欧美国产精品| 噜噜噜噜噜久久久久久91| 欧美高清成人免费视频www| 99在线人妻在线中文字幕| 午夜久久久久精精品| 日韩三级伦理在线观看| 久久这里有精品视频免费| 我要搜黄色片| 国产毛片a区久久久久| 午夜老司机福利剧场| 不卡一级毛片| 人妻系列 视频| 久久精品国产亚洲av香蕉五月| 国产午夜福利久久久久久| 天堂√8在线中文| 日日摸夜夜添夜夜添av毛片| 性插视频无遮挡在线免费观看| 身体一侧抽搐| 国产视频内射| av天堂在线播放| 日本在线视频免费播放|