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

    基于WRFDA-Chem 三維變分同化系統(tǒng)的中國PM2.5 和PM10 資料同化試驗(yàn)

    2022-10-21 03:24:48魏穎趙秀娟張自銀徐敬劉志權(quán)孫溦陳丹
    氣候與環(huán)境研究 2022年5期
    關(guān)鍵詞:氣溶膠偏差粒徑

    魏穎 趙秀娟 張自銀 徐敬 劉志權(quán) 孫溦 陳丹

    1 北京城市氣象研究院,北京 100089

    2 美國大氣研究中心,科羅拉多州博爾德 80301

    3 中國氣象科學(xué)研究院災(zāi)害天氣國家重點(diǎn)實(shí)驗(yàn)室,北京 100081

    1 引言

    近年來,中國的大城市面臨以細(xì)顆粒物(PM2.5)和高濃度臭氧(O3)為典型特征的大氣復(fù)合污染問題(朱彤等, 2010; Zhang and Cao, 2015;Sun and Chen, 2017; Zhao et al., 2020)。大氣顆粒物污染影響大氣能見度,并通過氣溶膠—輻射—云相互作用產(chǎn)生顯著的氣候效應(yīng)和生態(tài)環(huán)境效應(yīng)(Zhang et al., 2010; Boucher et al., 2014),同時(shí)影響人體健康(闞海東和陳秉衡, 2002; Liang et al.,2018)。當(dāng)前國際主流大氣化學(xué)模式對(duì)我國大氣污染模擬預(yù)報(bào)存在明顯偏差,秋冬季重污染期間PM2.5 濃度的預(yù)報(bào)偏差高達(dá)30%~50%(Zheng et al., 2015; Gao et al., 2016; Wu et al., 2020)。大氣化學(xué)模式的模擬結(jié)果與觀測(cè)存在較大偏差,其不確定性有多方面來源,從模式對(duì)污染形成和傳輸?shù)谋磉_(dá)角度,主要來自于以下3 個(gè)方面:1)排放清單的不確定性(Zheng et al., 2009; 曹國良等, 2011; 魏巍等, 2011);2)輸入氣象場(chǎng)的不確定性,尤其風(fēng)場(chǎng)、濕度、降水、溫度等氣象因子對(duì)大氣化學(xué)過程模擬會(huì)產(chǎn)生顯著影響(Zhang, 2008; Li et al., 2020);3)模式理化過程方案的不確定性(王穎等, 2010; 張穎等, 2016)。除了以上與污染形成和傳輸?shù)冗^程密切相關(guān)的3 個(gè)方面之外,數(shù)值模式技術(shù)本身的局限,如空間分辨率和理化方案適用尺度的匹配等,也會(huì)影響模擬結(jié)果的可靠性。減小空氣質(zhì)量預(yù)報(bào)不確定性對(duì)我國當(dāng)前大氣污染防控具有非常重要的現(xiàn)實(shí)意義和科學(xué)價(jià)值,在多方面不確定性來源混雜而無法快速診斷并提升模式性能的情況下,基于觀測(cè)約束的資料同化是提高預(yù)報(bào)效果的核心技術(shù)。

    早期限制大氣污染資料同化的主要因素是觀測(cè)資料的有限性(朱江等, 2018),近年隨著觀測(cè)數(shù)據(jù)增加,世界各國已開展了基于不同方法的大氣污染資料同化方法研究:1)較為簡(jiǎn)便易行的最優(yōu)插值法,如應(yīng)用于中國上海地區(qū)NO2、SO2及PM10的三維分析場(chǎng)預(yù)報(bào)(崔應(yīng)杰等, 2006)、歐洲地區(qū)PM10 的總濃度預(yù)報(bào)(Tombette et al., 2009)等;2)基于集合卡爾曼濾波系統(tǒng)的初始場(chǎng)、排放源或二者的同時(shí)同化,如基于嵌套網(wǎng)格空氣質(zhì)量預(yù)報(bào)模式系統(tǒng)(Nested Air Quality Prediction Modeling System,NAQPMS)開展的地面觀測(cè)對(duì)O3及其前體物的初始場(chǎng)和排放源同時(shí)同化(Tang et al., 2010),Miyazaki et al.(2012)同化NO2柱濃度資料得到全球NOx排放源,基于天氣大氣化學(xué)模式(Weather Research and Forecasting model coupled with Chemistry, WRF-Chem)開展對(duì)PM2.5 初始場(chǎng)和排放源的同時(shí)同化預(yù)報(bào)(Peng et al., 2017),基于同化衛(wèi)星資料的氣溶膠同化(Dai et al., 2019)等;3)基于全局優(yōu)化的三維變分同化技術(shù),如中國氣象科學(xué)研究院基于區(qū)域環(huán)境氣象數(shù)值模式(Chinese Unified Atmospheric Chemistry Environment,CUACE)開展的氣溶膠同化研究(Niu et al., 2008),Pagowski et al.(2010)對(duì)北美地區(qū)PM2.5 和O3濃度的同化預(yù)報(bào)試驗(yàn),Li et al.(2013)同化洛杉磯地面PM2.5 及其組分質(zhì)量濃度并進(jìn)行驗(yàn)證以及美國大氣研究中心開展的一系列基于格點(diǎn)統(tǒng)計(jì)插值(Gridpoint Statistical Interpolation, GSI)同化系統(tǒng)針對(duì)北美和中國地區(qū)的初始場(chǎng)同化工作(Liu et al., 2011; Jiang et al., 2013; Chen et al., 2019),汪代春等(2020)基于WRF-Chem 開展的中國氣溶膠和氣態(tài)前體物協(xié)同同化等;4)將空間全局最優(yōu)擴(kuò)展到空間和時(shí)間全局最優(yōu)的四維變分技術(shù),如Yumimoto et al.(2008)通過同化雷達(dá)消光系數(shù)反演起沙通量,Guerrette and Henze(2015)構(gòu)建伴隨模式實(shí)現(xiàn)對(duì)黑碳?xì)馊苣z的同化研究,Zhang et al.(2015)利用伴隨模式進(jìn)行了顆粒物污染源解析等。隨著研究的開展,同化技術(shù)也開始應(yīng)用于業(yè)務(wù)系統(tǒng)中,例如,中國環(huán)境監(jiān)測(cè)總站的NO2、SO2及PM2.5 空氣質(zhì)量預(yù)報(bào)采用最優(yōu)插值法同化,24 h 預(yù)報(bào)均方根誤差能夠減小23%(Zheng et al., 2018)。但總體來說在國內(nèi)關(guān)于大氣污染業(yè)務(wù)預(yù)報(bào)系統(tǒng)中的資料同化研究還較為缺乏。

    在資料同化方法中,三維變分將最小方差估計(jì)問題轉(zhuǎn)化為求代價(jià)函數(shù)最小值的問題(Sasaki,1970),可以考慮較為復(fù)雜和弱非線性的觀測(cè)算子,且計(jì)算代價(jià)小,往往作為常規(guī)業(yè)務(wù)資料同化的首選方法(Niu et al., 2008; Zhang et al., 2011)。由于氣溶膠多粒徑段、多化學(xué)成分組成的特點(diǎn),觀測(cè)資料形式多樣(如不同觀測(cè)手段、不同粒徑段和物種的氣溶膠觀測(cè)數(shù)據(jù)以及間接觀測(cè)產(chǎn)品氣溶膠光學(xué)厚度AOD),而不同模擬方案中的氣溶膠的粒徑/組分分類不同,因此面向不同的觀測(cè)系統(tǒng)和不同的大氣化學(xué)模式變量,必須建立不同的氣溶膠資料同化系統(tǒng)(Benedetti et al., 2009; Pagowski et al., 2010; Li et al., 2013; 靳璐濱等, 2016; Zang et al., 2016),增加了同化系統(tǒng)開發(fā)和移植的技術(shù)難度。

    為了應(yīng)對(duì)多組分多粒徑氣溶膠表達(dá)在模式中的復(fù)雜性,近年國內(nèi)外學(xué)者基于“一步法”(直接對(duì)所有的氣溶膠模式變量實(shí)施同時(shí)同化),在GSI同化平臺(tái)做了一系列工作,如Liu et al.(2011)針對(duì)GOCART 氣溶膠方案,在GSI 同化平臺(tái)上發(fā)展了對(duì)MODIS AOD 產(chǎn)品的三維變分同化算法,實(shí)現(xiàn)了對(duì)模式14 個(gè)化學(xué)變量的直接同化;Jiang et al.(2013)將WRF/Chem 中的GOCART 氣溶膠方案引入GSI 同化平臺(tái)進(jìn)行地面PM10 觀測(cè)資料同化;Schwartz et al.(2012)在同化MODIS 的AOD 基礎(chǔ)上新增了地面PM2.5 總量的三維變分同化;Chen et al.(2019)在GSI 同化平臺(tái)引入針對(duì)WRF/Chem中MOSAIC 氣溶膠機(jī)制的同化方案,實(shí)現(xiàn)對(duì)PM2.5和PM10 包含的32 個(gè)分粒徑段氣溶膠變量的直接同化,并將其應(yīng)用于中國冬季重污染成因分析。

    為了提升大氣污染預(yù)報(bào)和能見度預(yù)報(bào)的業(yè)務(wù)服務(wù)能力,并推動(dòng)大氣污染業(yè)務(wù)預(yù)報(bào)系統(tǒng)(RMAPSChem,WRF-Chem 內(nèi)核)與氣象業(yè)務(wù)預(yù)報(bào)系統(tǒng)(RMAPS-ST,WRF 內(nèi)核,氣象同化核心為WRFDA)的耦合,并減輕業(yè)務(wù)系統(tǒng)維護(hù)壓力,北京城市氣象研究院與美國大氣研究中心合作,在RMAPS-ST所采用的氣象同化平臺(tái)WRFDA 基礎(chǔ)上開發(fā)了與GSI 化學(xué)同化模塊類似的化學(xué)同化平臺(tái)(WRFDAChem),設(shè)計(jì)實(shí)現(xiàn)能靈活接入常用氣溶膠方案,提高氣溶膠初始場(chǎng)模擬準(zhǔn)確性,通過同化對(duì)比試驗(yàn)指明不同同化方案對(duì)氣溶膠模擬的影響,同時(shí)利用同化試驗(yàn)診斷發(fā)現(xiàn)模式對(duì)粗顆粒的顯著低估與排放信息表達(dá)缺失有關(guān)以及模式常用的非均相反應(yīng)攝取系數(shù)在新增同化前體物的情形下應(yīng)予以新的調(diào)整,以促進(jìn)模擬同化系統(tǒng)的協(xié)同發(fā)展(Sun et al., 2020)。本文基于RMAPS-Chem + WRFDA-Chem 同化預(yù)報(bào)系統(tǒng),以2016 年秋冬季為例進(jìn)行了為期1 個(gè)月的同化對(duì)照試驗(yàn),研究WRFDA-Chem 同化地面觀測(cè)資料對(duì)中國不同區(qū)域不同污染特點(diǎn)的PM2.5 和PM10預(yù)報(bào)的影響,同時(shí)對(duì)比分析了不同同化頻次對(duì)同化效果的影響;在此基礎(chǔ)上,結(jié)合模式本身的理化方案和偏差來源進(jìn)一步探討了在業(yè)務(wù)系統(tǒng)中開展同化預(yù)報(bào)的制約因素和發(fā)展方向,為大氣污染業(yè)務(wù)預(yù)報(bào)和同化的相關(guān)研究提供參考。

    2 模式與同化系統(tǒng)

    2.1 RMAPS-Chem 系統(tǒng)介紹

    睿圖—化學(xué)預(yù)報(bào)系統(tǒng)(RMAPS-Chem)是基于北京區(qū)域環(huán)境氣象數(shù)值預(yù)報(bào)系統(tǒng)(BREMPS)改進(jìn)發(fā)展而來的業(yè)務(wù)模式系統(tǒng),其采用全球通用的大氣化學(xué)模式WRF/Chem 作為主積分模式,即中尺度氣象模式(WRF)和大氣化學(xué)模式(Chem)在線耦合同步積分計(jì)算,并由采用快速更新同化技術(shù)的睿圖—短期預(yù)報(bào)系統(tǒng)(RMAPS-ST)提供氣象背景場(chǎng)。RMAPS-Chem 采用CBMZ 氣相化學(xué)機(jī)制和MOSAIC氣溶膠機(jī)制(未開啟高濕度下SO2-NO2-NO3通過快速非均相過程生成硫酸鹽—硝酸鹽—銨鹽的過程表達(dá))(Chen et al., 2016),經(jīng)過多次升級(jí),基于京津冀地區(qū)觀測(cè)數(shù)據(jù)進(jìn)行了能見度參數(shù)化方案改進(jìn)和SO2排放清單優(yōu)化(徐敬等, 2019; 趙秀娟等, 2019),近年來為北京重大活動(dòng)和華北地區(qū)重霾過程和污染氣象條件預(yù)報(bào)提供技術(shù)支撐。本研究所采用的RMAPSChem v2.0 是在RMAPS-Chem v1.0 的基礎(chǔ)上進(jìn)行配置更新和升級(jí)后的新版本,主要圍繞擴(kuò)大預(yù)報(bào)區(qū)域范圍、保證預(yù)報(bào)時(shí)效的目的,在水平網(wǎng)格范圍、垂直分層和模式版本等配置上進(jìn)行升級(jí),同時(shí)對(duì)排放清單和地表靜態(tài)數(shù)據(jù)集做了更新。RMAPSChem v2.0 目前使用WRF v4.0.2 版本,采用蘭勃特投影水平分辨率9/3 km 的嵌套網(wǎng)格,外層覆蓋我國東部大部分地區(qū)(即研究區(qū)域,見圖1),內(nèi)層覆蓋京津冀與汾渭平原地區(qū)。模式垂直分層為30 層,其中1500 m 以下為12 層,保障邊界層內(nèi)污染物排放、輸送擴(kuò)散和化學(xué)轉(zhuǎn)化等過程的精細(xì)化模擬和污染物的可靠預(yù)報(bào)。靜態(tài)排放源則引入了清華大學(xué)最新發(fā)布的0.1°分辨率的MEIC-2016 清單替代原來的業(yè)務(wù)清單(Zhang et al., 2019)。地表分類采用了基于2000~2010 年MODIS 遙感數(shù)據(jù)建立的Landuse 數(shù)據(jù)集,能更好地表征中國城市地表特征和陸地利用情況。在業(yè)務(wù)系統(tǒng)中,每天從00:00(協(xié)調(diào)世界時(shí),下同)循環(huán)熱啟以提供96 h預(yù)報(bào),初始場(chǎng)來自于前一天24 h 預(yù)報(bào)結(jié)果。

    圖1 研究區(qū)域范圍及同化試驗(yàn)中采用觀測(cè)站點(diǎn)(圓點(diǎn))位置。方框?yàn)? 個(gè)分析區(qū)域:華北平原(NCP)、西北地區(qū)(NWC)、長三角地區(qū)(YRD)Fig. 1 Model domains and locations of the observational sites (dots).The regions defined in rectangles represent the three regions used in the analysis: North China Plain (NCP), Northwest China (NWC), and Yangtze River Delta (YRD)

    2.2 WRFDA-Chem 同化系統(tǒng)

    北京城市氣象研究院聯(lián)合美國大氣研究中心(NCAR)在WRFDA 氣象同化平臺(tái)基礎(chǔ)上開發(fā)了全新的多物種多粒徑段氣溶膠3D-VAR 同化系統(tǒng)WRFDA-Chem(Sun et al., 2020),采用和現(xiàn)有氣象業(yè)務(wù)同化系統(tǒng)RMAPS-ST 完全相同的同化內(nèi)核(WRFDA),簡(jiǎn)化了未來耦合系統(tǒng)研發(fā)和業(yè)務(wù)運(yùn)行成本。在采用3D-VAR 同化初始場(chǎng)中,通過構(gòu)造代價(jià)函數(shù)J(x)將最小方差估計(jì)問題轉(zhuǎn)化為求代價(jià)函數(shù)最小值的問題,由此避免了直接求解權(quán)重矩陣,代價(jià)函數(shù)J(x)的表達(dá)式如下:

    其中,x為分析變量,xb為背景場(chǎng)矢量;B為背景誤差協(xié)方差矩陣,分析中決定背景場(chǎng)誤差和觀測(cè)誤差權(quán)重分配;y為觀測(cè)場(chǎng)矢量;R是觀測(cè)誤差協(xié)方差矩陣,H為觀測(cè)算子。其中,對(duì)背景場(chǎng)誤差協(xié)方差矩陣B的估計(jì)問題是資料同化的核心,它決定了同化時(shí)觀測(cè)信息對(duì)分析變量(即各個(gè)組分/粒徑變量)的訂正范圍和程度。WRFDA-Chem 同化系統(tǒng)通過水平、垂直和物理三方面的控制變量變換實(shí)現(xiàn)簡(jiǎn)化處理,并通過標(biāo)準(zhǔn)方差和水平/垂直相關(guān)尺度參數(shù)反映各變量之間的具體相關(guān)和結(jié)構(gòu);為B矩陣服務(wù)的氣溶膠和氣體物種變量的相關(guān)統(tǒng)計(jì)參數(shù)(Descombes et al., 2015),通過GEN_BE v2.0 模塊計(jì)算,采用NMC 方法(Parrish and Derber, 1992)生成。本文基于研究時(shí)段2016 年11 月為期1 個(gè)月的WRF/Chem 預(yù)報(bào),計(jì)算B統(tǒng)計(jì)量。

    針對(duì)顆粒物的多種觀測(cè)資料,WRFDA-Chem系統(tǒng)設(shè)計(jì)了多種同化方案以表達(dá)不同的觀測(cè)算子,既可單獨(dú)同化PM2.5 和PM10,也可同時(shí)同化二者;在同時(shí)同化的過程中,為了更好地表達(dá)粒徑間關(guān)系,還設(shè)計(jì)了針對(duì)PM2.5 和PMcoarse(粗顆粒,PM10-2.5)觀測(cè)算子的同時(shí)同化方案。在顆粒物觀測(cè)資料外,WRFDA-Chem 中還加入了針對(duì)4 種常規(guī)氣態(tài)污染物(SO2、NO2、CO、O3)觀測(cè)資料的觀測(cè)算子,能實(shí)現(xiàn)對(duì)顆粒物觀測(cè)和氣態(tài)污染物觀測(cè)的同時(shí)同化,以更好地表達(dá)氣體污染物對(duì)二次生成顆粒物的過程。為了適應(yīng)不同氣溶膠方案,WRFDAChem 設(shè)計(jì)采用注冊(cè)表輸入方式靈活切換分析變量;目前已接入了GOCART 和MOSAIC 兩種氣溶膠機(jī)制。

    與RMAPS-Chem 匹配,本研究中同化過程所采用的氣溶膠分析變量是適用于中國地區(qū)尤其是人為源大氣污染狀況的MOSAIC 機(jī)制,氣溶膠成分包含硫酸鹽(SO)、硝酸鹽(N)、銨鹽(N)、氯鹽(NA)、鈉鹽(Cl)、黑碳(BC)、有機(jī)碳(OC)和其他無機(jī)成分(OIN)等8 個(gè)物種,每個(gè)物種采用0.039~0.1、0.1~1.0、1.0~2.5 和2.5~10 μm 4 個(gè)分段粒徑,共計(jì)32 個(gè)變量。其中顆粒物同化過程中,PM2.5 總質(zhì)量濃度同化影響前3 個(gè)粒徑段8 種氣溶膠物種,共計(jì)24 個(gè)變量(見公式2),PM10總質(zhì)量濃度觀測(cè)的影響則是通過PMcoarse 的觀測(cè)算子,影響第4 個(gè)粒徑段共計(jì)8 個(gè)變量(見公式3和4)。同化時(shí)對(duì)觀測(cè)PM2.5質(zhì)量濃度的分配采用直接同化的“一步法”,在考慮背景誤差和觀測(cè)誤差的情況下,通過對(duì)背景誤差協(xié)方差的統(tǒng)計(jì),考慮多物種多粒徑段組分為變量的整體代價(jià)函數(shù)最小化計(jì)算得來:

    其中ρd為干空氣密度。

    在計(jì)算代價(jià)函數(shù)最小值時(shí),對(duì)觀測(cè)誤差協(xié)方差的估計(jì)考慮儀器誤差和觀測(cè)代表性誤差,其中觀測(cè)代表性誤差估計(jì)采用下式計(jì)算:

    其中,γ 為比例因子,設(shè)定為0.5;ε0為儀器誤差,為1.0+0.0075Mi,Mi為污染物濃度(單位:μg m-3);Δx為網(wǎng)格分辨率,本文為9 km;L為觀測(cè)影響半徑,設(shè)定為2 km。上述參數(shù)根據(jù)Chen et al.(2019)的同化試驗(yàn)敏感性測(cè)試結(jié)果設(shè)定。

    在利用地面觀測(cè)數(shù)據(jù)進(jìn)行大氣污染同化時(shí),單個(gè)觀測(cè)點(diǎn)與模式網(wǎng)格面不匹配是在同化結(jié)果中引入誤差的重要因素。當(dāng)前對(duì)觀測(cè)代表性誤差的估計(jì)一般考慮網(wǎng)格分辨率、觀測(cè)站點(diǎn)類別和物種等要素的影響(Elbern et al., 2007)。有研究利用京津冀地面觀測(cè)網(wǎng)絡(luò)對(duì)大氣污染物地面觀測(cè)代表誤差進(jìn)行了估計(jì),結(jié)果表明網(wǎng)格內(nèi)PM2.5 和PM10 的代表性誤差分別為11.1 μg m-3和22.2 μg m-3(不確定性11%)(李飛等, 2019)。本文暫未針對(duì)觀測(cè)代表性誤差對(duì)同化結(jié)果的影響進(jìn)行單獨(dú)評(píng)估,將在未來的工作中考慮進(jìn)行相關(guān)敏感性試驗(yàn)評(píng)估其影響。

    3 背景誤差協(xié)方差統(tǒng)計(jì)分析

    氣溶膠同化的難點(diǎn)就是利用有限散點(diǎn)的PM2.5地面總質(zhì)量濃度資料對(duì)模式中多物種多粒徑段氣溶膠組分進(jìn)行約束。如前文所述,是通過控制變量變換簡(jiǎn)化處理,通過水平、垂直和物理三方面的轉(zhuǎn)化實(shí)現(xiàn)簡(jiǎn)化,通過標(biāo)準(zhǔn)方差和水平/垂直相關(guān)尺度參數(shù)反映各變量之間的具體相關(guān)和結(jié)構(gòu)。我們采用了NMC 方法生成了背景誤差協(xié)方差矩陣中簡(jiǎn)化的各統(tǒng)計(jì)參數(shù)。圖2 為全區(qū)域平均的各變量的水平相關(guān)尺度參數(shù),可以看到4 個(gè)粒徑段內(nèi)的氣溶膠各組分,除量級(jí)較低的鈉鹽(Na+)和氯鹽(Cl-)外,其他組分(硫酸鹽SO、硝酸鹽N、銨鹽N、黑碳BC、有機(jī)碳OC、其他無機(jī)成分OIN)的水平相關(guān)尺度在地面附近大多在15~25 km 左右,結(jié)合模式網(wǎng)格9 km 的分辨率,以上水平相關(guān)尺度較為符合預(yù)期。在對(duì)流層內(nèi),各變量的水平相關(guān)尺度隨著高度增加逐步減小。圖3 給出了全區(qū)域平均的各變量的模式背景誤差標(biāo)準(zhǔn)偏差,可以看到不同粒徑段內(nèi)不同組分的背景誤差的標(biāo)準(zhǔn)偏差存在很大的差異。在構(gòu)成PM2.5 的小粒徑段內(nèi),硝酸鹽具有最大的標(biāo)準(zhǔn)偏差,其次是銨鹽,地面附近OIN 也具有較大的標(biāo)準(zhǔn)偏差,這反映了模式在沒考慮新增非均相化學(xué)過程時(shí),模擬區(qū)域中PM2.5 濃度的表達(dá)以硝酸鹽成分為主的情況,而在大粒徑段內(nèi)OIN 組分的標(biāo)準(zhǔn)偏差量級(jí)顯著大于其他組分,這是由于在大粒徑段氣溶膠成分以沙塵等自然源為主,二次無機(jī)鹽組分占比較小。由于大部分排放源分布在地面層,氣溶膠質(zhì)量濃度在邊界層以上逐漸降低,因而各粒徑段各物種變量的背景誤差標(biāo)準(zhǔn)差隨著垂直高度增加而降低,在對(duì)流層頂部下降至接近0。上述氣溶膠組分在不同粒徑段的標(biāo)準(zhǔn)差的分布變化體現(xiàn)了中國中東部地區(qū)顆粒物由多組分多粒徑段氣溶膠構(gòu)成的特點(diǎn)。

    圖2 2016 年11 月模式區(qū)域顆粒物粒徑段(a)bin1、(b)bin2、(c)bin3 和(d)bin4 控制變量的水平相關(guān)尺度(單位:km)Fig. 2 Length scale of the aerosol species (units: km) of the (a) first size bin, (b) second size bin, (c) third size bin, and (d) fourth size bin over the domain in Nov 2016

    圖3 同圖2,但為PM2.5 控制變量的背景誤差標(biāo)準(zhǔn)差垂直廓線(單位:μg m-3)Fig. 3 Background error standard deviations (μg m-3) of the aerosol species of the (a) first size bin, (b) second size bin, (c) third size bin, and (d)fourth size bin over the domain in Nov 2016

    4 資料與方法

    4.1 模式設(shè)置與試驗(yàn)設(shè)計(jì)

    本研究針對(duì)RMAPS-Chem 外層9 km×9 km區(qū)域進(jìn)行同化試驗(yàn),范圍為(24°N~46°N,97°E ~122°E),覆蓋中國大陸大部分區(qū)域,網(wǎng)格數(shù)為304×304(研究區(qū)域范圍詳見圖1)。研究時(shí)間段為2016 年11 月,其中初始場(chǎng)模擬的預(yù)熱時(shí)間為2016 年10 月24~31 日共計(jì)8 d。氣象場(chǎng)初始條件和邊界條件采用美國環(huán)境預(yù)報(bào)中心(NCEP)的1°(緯度)×1°(經(jīng)度)FNL 全球再分析數(shù)據(jù),化學(xué)場(chǎng)初始條件和邊界條件采用WRF-Chem 提供的默認(rèn)方案。模式采用的主要物理方案信息見表1,人為排放源采用清華大學(xué)0.1°(緯度)×0.1°(經(jīng)度)的MEIC-2016 清單(Zhang et al., 2019)經(jīng)插值匹配到RMAPS-Chem 的分辨率,包含甲烷、一氧化碳、二氧化硫(SO2)、氮氧化物(NOx)、揮發(fā)性有機(jī)物、氨、一次顆粒物(含元素碳、有機(jī)碳)等30 個(gè)物種。由于人為源相關(guān)的無組織揚(yáng)塵排放不僅受自然條件影響,同時(shí)也受人為活動(dòng)和減排管控措施影響較大,準(zhǔn)確的排放總量統(tǒng)計(jì)和時(shí)空分配仍然存在較大困難,為了保證業(yè)務(wù)運(yùn)行的一致性和穩(wěn)定性,在業(yè)務(wù)模式中尚沒有考慮人為源無組織揚(yáng)塵排放。另外,由于現(xiàn)有模式的自然源沙塵模塊中,采用了半經(jīng)驗(yàn)的起沙參數(shù)化方案,起沙過程并未對(duì)不同起沙機(jī)制(空氣動(dòng)力直接拖曳作用、沙粒的躍移撞擊、土壤聚合體的分裂)進(jìn)行區(qū)分(張宏升和李曉嵐, 2014),也未考慮不同粒徑的沙塵粒子微觀運(yùn)動(dòng)方式(蠕移、躍移、懸浮運(yùn)動(dòng))的差異,對(duì)于起沙參數(shù)(如臨界起沙摩擦速度、土壤濕度等)并沒有開展很好的本地化驗(yàn)證,對(duì)垂直方向的抬升過程和分布表達(dá)也較為簡(jiǎn)易,容易造成虛假起沙和局地近地面層沙塵通量過大,從而帶來業(yè)務(wù)模擬效果的不穩(wěn)定,所以在業(yè)務(wù)系統(tǒng)中也未開啟WRF-Chem 模式的沙塵自然源排放過程。

    為了分析同化地面觀測(cè)資料對(duì)顆粒物預(yù)報(bào)的影響,本研究設(shè)置了同化和不同化2 組對(duì)比試驗(yàn),同化試驗(yàn)采用逐6 小時(shí)循環(huán)同化(DA_6h),另一組則采用不同化(NODA)作為對(duì)照試驗(yàn)。同時(shí)為了探究不同頻次的同化對(duì)預(yù)報(bào)的影響,另外設(shè)計(jì)了逐3 小時(shí)循環(huán)同化(DA_3h)和逐1 小時(shí)循環(huán)同化(DA_1h)2 組試驗(yàn)。研究采用同時(shí)同化顆粒物PM2.5 和PMcoarse 的方案,并用同化得到的00:00的結(jié)果作為初始場(chǎng)進(jìn)行24 h 預(yù)報(bào)試驗(yàn),探討同化初始場(chǎng)對(duì)于污染預(yù)報(bào)的作用,NODA 則采用原來預(yù)報(bào)循環(huán)熱啟的結(jié)果作為初始場(chǎng),同化方案和試驗(yàn)設(shè)置詳見表1。

    表1 模擬采用的參數(shù)化方案和試驗(yàn)設(shè)計(jì)Table 1 Information regarding the parameterization and the detailed settings of the four experiments

    4.2 觀測(cè)資料與統(tǒng)計(jì)

    本文觀測(cè)數(shù)據(jù)來自中國環(huán)境監(jiān)測(cè)總站發(fā)布的國控站常規(guī)6 項(xiàng)污染物濃度觀測(cè)結(jié)果,覆蓋研究區(qū)域共計(jì)836 個(gè)站點(diǎn)(分布見圖3),觀測(cè)時(shí)間段為2016 年11 月。為了保證觀測(cè)資料的質(zhì)量,同化前首先對(duì)觀測(cè)數(shù)據(jù)進(jìn)行了質(zhì)量控制,包括極值控制和異常值剔除:對(duì)PM2.5 大于500 μg m-3,PM10 大于800 μg m-3的數(shù)據(jù)進(jìn)行剔除;同時(shí)剔除結(jié)果小于0 以及PM10 觀測(cè)小于PM2.5 的異常值。為保證同化的穩(wěn)定運(yùn)行,在同化過程中對(duì)與背景值(模式模擬值投射到觀測(cè)空間)相比差異大于120 μg m-3的觀測(cè)值也進(jìn)行了剔除。

    在定量化評(píng)估WRFDA-Chem 同化系統(tǒng)對(duì)預(yù)報(bào)的效果改善情況時(shí),引入平均偏差(MB),歸一化平均偏差(NMB),均方根誤差(RMSE)和相關(guān)系數(shù)(R)等統(tǒng)計(jì)指標(biāo)進(jìn)行檢驗(yàn)分析,且僅在觀測(cè)有效數(shù)據(jù)大于總時(shí)次的2/3 時(shí)認(rèn)為該站點(diǎn)為有效站點(diǎn)。統(tǒng)計(jì)時(shí)根據(jù)觀測(cè)站點(diǎn)經(jīng)緯度插值提取相近格點(diǎn)的模擬值,將區(qū)域內(nèi)所有有效的單個(gè)站點(diǎn)的觀測(cè)值和模擬值按時(shí)間序列依次接續(xù)作為樣本計(jì)算相關(guān)的統(tǒng)計(jì)參數(shù),作為區(qū)域平均值。

    5 結(jié)果分析

    5.1 分析場(chǎng)效果

    對(duì)大氣化學(xué)數(shù)值模式而言,初始場(chǎng)的準(zhǔn)確性直接影響模式預(yù)報(bào)的效果,下面對(duì)同化前/后00:00 的初始場(chǎng)/分析場(chǎng)進(jìn)行對(duì)比分析。圖4a 和4b 分別給出了模式同化前的00:00 初始場(chǎng)(來自于前一天00:00 熱啟的24 h 預(yù)報(bào)結(jié)果)PM2.5 和PM10 月平均值以及與觀測(cè)值的對(duì)比。如圖所示,模式對(duì)我國北方大部分區(qū)域尤其是華北的PM2.5 濃度整體預(yù)報(bào)能力較好,京津冀—山東—河南污染帶和汾渭平原地區(qū)的PM2.5 濃度在100~200 μg m-3之間,河西地區(qū)和京津冀北部城市PM2.5 濃度在30~100 μg m-3之間。然而模式對(duì)我國中南部地區(qū)的PM2.5 濃度存在明顯的高估,例如對(duì)湖南、湖北地區(qū)高估超過50 μg m-3,對(duì)四川盆地的高估則超過100 μg m-3,可能反映了排放清單對(duì)這些區(qū)域的冬季排放量存在高估,或者模式本身對(duì)四川盆地的復(fù)雜地形和冬季污染表達(dá)仍然有較大不確定性。PM10方面,由于我國北方下墊面(地表粗糙度、可風(fēng)蝕程度、土壤濕度等)和空氣動(dòng)力學(xué)(風(fēng)速)條件比南方更利于起沙,而且降水相對(duì)較少,除了春季的長距離傳輸?shù)淖匀簧硥m天氣過程,也易于產(chǎn)生小范圍起沙和局地?zé)o組織揚(yáng)塵,即使在秋冬季節(jié)我國北方觀測(cè)的PM10 濃度也達(dá)到了150 μg m-3以上。而這部分小范圍起沙和局地?zé)o組織揚(yáng)塵,不僅受自然條件影響,同時(shí)也受人為活動(dòng)和減排管控措施影響,其排放總量和時(shí)空分配都存在較大不確定性;目前在模式的靜態(tài)排放源和起沙模塊中都沒有表達(dá),因此模式對(duì)PM10 濃度的北方特異性很難準(zhǔn)確把握。如圖4d 所示,大部分北方地區(qū)模擬低估0~100%(約100~200 μg m-3),而中南部地區(qū)則高估20%~100%(約100~200 μg m-3)。由于人為源排放、自然源過程、以及模式本身理化過程的多方面不確定性,目前模式對(duì)于南北方的PM10 偏差存在方向性差異。

    圖4 2016 年11 月模式區(qū)域分析時(shí)刻(00:00)未同化的(a)PM2.5 和(b)PM10 月平均濃度(陰影)與地面站點(diǎn)(圓點(diǎn))觀測(cè)值對(duì)比、同化前00:00(c)PM2.5 和(d)PM10 月平均濃度(陰影)歸一化平均偏差、同化前后(DA_6h-NODA)(e)PM2.5 和(f)PM10 月平均濃度(陰影)差異Fig. 4 Surface distributions of monthly mean (a) PM2.5 and (b) PM10 concentrations (shandings) at 0000UTC for NODA against the observations(circles), the normalized mean bias of (c) PM2.5 and (d) PM10 concentrations (shandings) at 0000UTC for NODA, and the difference of (e) PM2.5 and(f) PM10 monthly mean concentrations (shandings) between DA_6h and NODA over the domain in Nov 2016

    經(jīng)過DA_6h 同化得到的PM2.5 濃度00:00 分析場(chǎng)如圖4e 所示,同化后的PM2.5 濃度在空間分布上得到調(diào)整,西北部分城市PM2.5 濃度提升20~50 μg m-3,中國南方大部分地區(qū)降低20~100 μg m-3,結(jié)合圖4a 可見同化后的PM2.5 濃度更接近觀測(cè)結(jié)果。圖4f 表明,同化系統(tǒng)對(duì)PM10模擬的改善作用也同樣明顯,同化地面觀測(cè)信息對(duì)我國北方地區(qū)PM10 濃度的提升可超過100 μg m-3,對(duì)南方PM10 濃度高估的問題也能降低50 μg m-3左右。圖5 和表2 分別給出了同化前后分析時(shí)刻的觀測(cè)—模擬散點(diǎn)對(duì)比和統(tǒng)計(jì)參數(shù)。如圖5a 所示,同化后模式對(duì)不同站點(diǎn)的高估和低估現(xiàn)象都得到明顯改善,DA_6h 將初始場(chǎng)的模擬偏差減少超過50%,NMB 和RMSE 也均減小一半左右。與此同時(shí),分析時(shí)刻的PM2.5 濃度相關(guān)性也增加了0.27,說明WRFDA-Chem 同化系統(tǒng)對(duì)污染物在分析時(shí)刻的時(shí)空相關(guān)性能起到明顯的調(diào)整作用。從圖5b 的散點(diǎn)圖來看,PM10 濃度的空間差異更大,同化后能夠減小這種差異,尤其是對(duì)PM10 濃度低估的問題有明顯改善。雖然由于區(qū)域內(nèi)東西部區(qū)域的偏差方向不同導(dǎo)致同化前后區(qū)域平均偏差變化較小,但是同化過程能顯著降低PM10 濃度均方根誤差30%~40%左右,且相關(guān)性顯著提升0.37。

    表2 同化前后各試驗(yàn)分析時(shí)刻PM2.5 和PM10 濃度的模擬統(tǒng)計(jì)參數(shù)Table 2 Statistics between model results and observations for PM2.5 and PM10 for different experiments at 0000UTC

    圖5 2016 年11 月模式區(qū)域同化前后分析時(shí)刻(00:00)的模擬與觀測(cè)(a)PM2.5 和(b)PM10 月平均濃度散點(diǎn)圖Fig. 5 Scatter plots between the model results and the observations for monthly mean concentration of (a) PM2.5 and (b) PM10 at 0000UTC for the experiments of NODA and DA_6h over the domain in Nov 2016

    需要指出的是,考慮到業(yè)務(wù)同化運(yùn)行的穩(wěn)定性,本研究中采用同化系統(tǒng)內(nèi)外兩重質(zhì)控來篩選觀測(cè)數(shù)據(jù),質(zhì)控過程的不同可能影響同化效果。第一重質(zhì)控是在進(jìn)同化系統(tǒng)前,對(duì)PM2.5 和PM10 觀測(cè)分別大于500 和800 μg m-3的觀測(cè)進(jìn)行剔除;第二重質(zhì)控是在同化過程中,對(duì)偏離背景值(模式模擬值投射到觀測(cè)空間)超過120 μg m-3的觀測(cè)進(jìn)行了剔除。由于某些地區(qū)(如京津冀地區(qū)、山西省、汾渭平原)的PM2.5 模擬結(jié)果在污染過程嚴(yán)重偏低,帶來部分觀測(cè)數(shù)據(jù)在第二重質(zhì)控過程中被剔除而無法用于同化過程,導(dǎo)致PM2.5 同化效果在這些地區(qū)不明顯。而由于PM10 觀測(cè)數(shù)據(jù)的同化是通過PMcoarse 影響第四個(gè)粒徑段的變量來開展;對(duì)PMcoarse 而言,120 μg m-3的剔除標(biāo)準(zhǔn)已較為寬泛,剔除數(shù)據(jù)量較少,因而PM10 在京津冀地區(qū)的同化效果優(yōu)于PM2.5。后期需要對(duì)雙重質(zhì)控,特別是第二重質(zhì)控中的剔除標(biāo)準(zhǔn)開展進(jìn)一步測(cè)試,以進(jìn)一步提高PM2.5 的同化效果。

    5.2 24 小時(shí)預(yù)報(bào)效果

    模式預(yù)報(bào)的準(zhǔn)確性不僅受到初始場(chǎng)的影響,同時(shí)也受排放源和氣象場(chǎng)的不確定性影響,除此之外,模式的系統(tǒng)性誤差也是影響預(yù)報(bào)結(jié)果的重要因素,因此隨著預(yù)報(bào)時(shí)效的增加,同化初始場(chǎng)對(duì)于改善預(yù)報(bào)的正向作用會(huì)逐漸減弱。為了檢驗(yàn)WRFDAChem 同化初始場(chǎng)對(duì)預(yù)報(bào)的改進(jìn)效果和持續(xù)影響時(shí)間,采用6 h 循環(huán)同化前后獲得的00:00 初始場(chǎng)進(jìn)行了為期1 個(gè)月的逐日預(yù)報(bào),并對(duì)0~24 h 預(yù)報(bào)的平均結(jié)果進(jìn)行了對(duì)比。如圖6 所示,24 h 以內(nèi)WRFDAChem 同化系統(tǒng)對(duì)PM2.5 濃度預(yù)報(bào)偏差的調(diào)整能降低10~25 μg m-3左右,RMSE 能降低20~40 μg m-3,相關(guān)性提高~0.3。表3 給出了24 h 內(nèi)平均統(tǒng)計(jì)參數(shù),同化后預(yù)報(bào)偏差降低50%以上,RMSE 降低25%,相關(guān)性提升14%。相較PM2.5 濃度而言,同化對(duì)PM10 濃度預(yù)報(bào)偏差的改善相對(duì)弱一些。由于PM10 濃度模擬值在我國北方和南方的偏差方向相反,同化后整個(gè)研究區(qū)域的平均預(yù)報(bào)偏差反而有所增加,但同化后24 h 內(nèi)平均RMSE 降低10%以上,且相關(guān)性明顯提升,超過25%??傮w而言,同化對(duì)PM2.5 和PM10 濃度預(yù)報(bào)效果的改善能持續(xù)24 h 以上,前12 小時(shí)改善效果更為明顯,其中PM2.5 的RMSE 能降低30%~50%。

    圖6 同化前、后的整個(gè)模擬區(qū)域(a1-a3)PM2.5 濃度和(b1-b3)PM10 濃度0~24 h 預(yù)報(bào)統(tǒng)計(jì)參數(shù)比較Fig. 6 Averaged bias, RMSE, and correlation for the simulated concentrations of (a1-a3) PM2.5 and (b1-b3) PM10 against the surface observations in different experiments over the 0-24 forecast hours

    表3 同化前后各試驗(yàn)PM2.5 和PM10 在0~24 h 預(yù)報(bào)中的統(tǒng)計(jì)參數(shù)Table 3 Statistics over the 0-24 forecast hours between model results and observations for PM2.5 and PM10 for different experiments

    由于顆粒物濃度的空間分布差異較大,且觀測(cè)站點(diǎn)東密西疏、空間分布不均,在同化系統(tǒng)迭代求解代價(jià)函數(shù)最小值進(jìn)行全局最優(yōu)化時(shí),對(duì)不同地區(qū)的同化效果也可能存在較大的差異。如圖7 所示,同化后華北平原、華中和四川盆地的PM2.5 濃度均方差顯著降低約30~60 μg m-3,華中地區(qū)和預(yù)報(bào)偏差較大的西北地區(qū)PM10 濃度均方差降低約20~40 μg m-3。為進(jìn)一步分析這一影響,本文對(duì)研究區(qū)域內(nèi)較受關(guān)注的華北平原、長三角地區(qū)和西北地區(qū)的24 h 預(yù)報(bào)進(jìn)行了統(tǒng)計(jì)分析(圖8)。未同化時(shí),華北平原PM2.5 濃度模擬偏差隨著預(yù)報(bào)時(shí)效(0~24 h)增加經(jīng)歷了高估—低估—高估的變化,這反映了在00:00 至24:00 的日內(nèi)變化過程中,原模式模擬在排放、氣象場(chǎng)、化學(xué)過程等多方面產(chǎn)生偏差的過程中,可能出現(xiàn)了正/負(fù)兩個(gè)不同方向的影響在不同時(shí)次存在協(xié)同和制約兩種相互作用的現(xiàn)象。由于初始場(chǎng)偏差和中間時(shí)次預(yù)報(bào)場(chǎng)偏差相反,導(dǎo)致同化后的預(yù)報(bào)偏差在某些時(shí)次反而有所增大。而長三角和西北地區(qū)的PM2.5 濃度偏差方向單一,能看到圖8b 和8c 中初始時(shí)刻的同化對(duì)減小0~24 h預(yù)報(bào)整體偏差能起到不同程度的正向作用,其中長三角地區(qū)的偏差改善更為明顯,24 h 預(yù)報(bào)時(shí)效內(nèi)偏差減小20%~50%,RMSE 減小12%~80%。對(duì)PM10而言,西北地區(qū)是重要的粗顆粒物源地,同化對(duì)于調(diào)整該地的PM10 濃度低估具有更為明顯的效果。相較長三角地區(qū)來說,同化系統(tǒng)對(duì)我國北方的PM10濃度預(yù)報(bào)改善效果更好,偏差減小更明顯,相關(guān)性提升的優(yōu)勢(shì)也能保持超過12 h。

    圖7 同化前(左)、DA_6h 同化后(右)24 h 預(yù)報(bào)與地面站點(diǎn)觀測(cè)(a、b)PM2.5 和(c、d)PM10 濃度的均方根誤差對(duì)比Fig. 7 RMSE maps at the surface validation sites for the concentration of (a, b) PM2.5 and (c, d) PM10 for 24-h forecasts of NODA (left column) and DA_6h (right column)

    圖8 同化前、后華北平原(左列)、長三角地區(qū)(中列)、西北地區(qū)(右列)(a1-a3、b1-b3、c1-c3)PM2.5 濃度和(d1-d3、e1-e3、f1-f3)PM10 濃度0~24 h 預(yù)報(bào)檢驗(yàn)統(tǒng)計(jì)參數(shù)Fig. 8 Averaged bias, RMSEs, and correlation coefficients for the simulated concentrations of (a1-a3, b1-b3, c1-c3) PM2.5 and (d1-d3, e1-e3,f1-f3) PM10 in the North China Plain (left colunm), the Yangtze River Delta (middle panel), and the Northwest China (right panel) against the surface observations in different experiments over the 0-24 forecast hours

    需要指出的是,在利用觀測(cè)同化來更新初始場(chǎng)的過程中,其本身的物理意義是將模式在一段時(shí)間內(nèi)來自多方面不確定性產(chǎn)生的累積偏差,在這一時(shí)刻予以強(qiáng)迫修正;而由于模式多方面不確定性(排放、氣象場(chǎng)、化學(xué)過程)本身在預(yù)報(bào)過程中(這里體現(xiàn)為日內(nèi)變化上)可能存在正/負(fù)兩個(gè)不同方向的協(xié)同和制約作用,因而在某一時(shí)刻對(duì)濃度場(chǎng)的強(qiáng)迫修正,可能帶來整體累積偏差在后續(xù)預(yù)報(bào)中反而增大的現(xiàn)象;在這種情況下,選擇哪個(gè)時(shí)次對(duì)濃度場(chǎng)進(jìn)行修正可能影響后續(xù)預(yù)報(bào)效果。整體而言,同化預(yù)報(bào)效果的提升其實(shí)與模式本身性能息息相關(guān),需要在深入理解模式不確定性和偏差來源的情況下,開展模式和同化系統(tǒng)的協(xié)同發(fā)展。由于本研究主要目的是服務(wù)于業(yè)務(wù)系統(tǒng),因而在試驗(yàn)設(shè)計(jì)上僅考慮了采用00:00 的同化分析場(chǎng)作為初始場(chǎng)來開展后續(xù)24 h 預(yù)報(bào),而未考慮其他時(shí)次,未來需要開展更多的試驗(yàn)以優(yōu)化同化預(yù)報(bào)流程。

    5.3 同化頻次對(duì)同化效果的影響

    對(duì)3D-VAR 而言,同化頻次也是影響同化效果的重要因素。為了探究RMAPS-Chem + WRFDAChem 同化預(yù)報(bào)的最佳方案,將3 h 和1 h 循環(huán)同化(DA_3h、DA_1h)與6 h 循環(huán)同化(DA_6h)的00:00 分析場(chǎng)及24 h 預(yù)報(bào)場(chǎng)進(jìn)行了對(duì)比分析(見圖6)。3 套循環(huán)同化試驗(yàn)在00:00 同化分析時(shí)刻的統(tǒng)計(jì)對(duì)比見表2,基于00:00 的初始場(chǎng)/同化分析場(chǎng)開展的24 h 預(yù)報(bào)結(jié)果見圖6 和表3。對(duì)比表明,更高的同化頻次能夠通過逐步減小背景場(chǎng)與觀測(cè)的差異從而進(jìn)一步降低00:00 時(shí)刻PM2.5 分析場(chǎng)的偏差,但是PM10 分析場(chǎng)偏差變化不大,DA_1h 中PM2.5 和PM10 的相關(guān)性都得到進(jìn)一步提升。隨著同化頻次提高對(duì)分析場(chǎng)的改進(jìn),PM2.5 的24 h 預(yù)報(bào)改進(jìn)也較為明顯,DA_1h 的24 h 預(yù)報(bào)偏差(MB、NMB)比DA_6h 減小了將近50%,相關(guān)性略提升0.01。PM10 方面,基于更高頻次的同化初始場(chǎng)開展的24 h 預(yù)報(bào)相關(guān)性有所提高,DA_1h 相關(guān)性較DA_6h 提高0.02,但偏差變化不大,這是由于模式對(duì)PM10 模擬在我國南北方偏差方向明顯相反,對(duì)整個(gè)區(qū)域模擬偏差的統(tǒng)計(jì)上可能存在抵消。而如圖6 所示,相較于DA_6h,DA_1h 能夠進(jìn)一步降低PM2.5 和PM10 的均方差,對(duì)12 h 后的預(yù)報(bào)相關(guān)性也能進(jìn)一步提高。因此在業(yè)務(wù)同化的應(yīng)用中,采用逐小時(shí)循環(huán)同化高頻次地引入觀測(cè)信息,可以獲得更為準(zhǔn)確的初始場(chǎng),從而有利于更長時(shí)效的提高預(yù)報(bào)準(zhǔn)確性。

    6 總結(jié)與討論

    為了服務(wù)于未來的氣象—化學(xué)一體化業(yè)務(wù)系統(tǒng),北京城市氣象研究院與NCAR 合作,在WRFDA氣象同化平臺(tái)基礎(chǔ)上新增了3D-VAR 大氣化學(xué)三維變分同化系統(tǒng)WRFDA-Chem,實(shí)現(xiàn)了對(duì)多物種多粒徑段氣溶膠觀測(cè)的同化功能。本研究將WRFDAChem 引入RMAPS-Chem 環(huán)境氣象數(shù)值預(yù)報(bào)系統(tǒng),利用地面顆粒物濃度觀測(cè)資料開展了多組同化試驗(yàn)的對(duì)比分析研究,結(jié)論如下:

    (1)利用WRFDA-Chem 同化地面PM2.5 和PM10 總質(zhì)量濃度觀測(cè),對(duì)污染物分析場(chǎng)的時(shí)空相關(guān)性能起到明顯的調(diào)整作用,DA_6h 同化后PM2.5和PM10 濃度相關(guān)性分別提升0.27 和0.37,同時(shí)有效減小污染物的空間分布偏差,同化后PM2.5濃度的平均偏差(MB)和歸一化平均偏差(NMB)降低50%以上,RMSE(均方根誤差)減小40%,對(duì)PM10 濃度的RMSE 也減少約40%。

    (2)WRFDA-Chem 通過同化初始濃度場(chǎng)對(duì)顆粒物預(yù)報(bào)效果的改進(jìn)能持續(xù)24 h 以上,PM2.5預(yù)報(bào)偏差降低50%以上,RMSE 降低25%,相關(guān)性提升14%,PM10 同化后24 h 內(nèi)平均RMSE 降低10%以上,相關(guān)性提升超過25%。同化在前12小時(shí)改善效果更為明顯,PM2.5 的RMSE 能降低30%~50%。

    (3)在具有不同污染特征的區(qū)域,同化效果存在較為明顯的差異,同化地面PM2.5 觀測(cè)對(duì)長三角地區(qū)的PM2.5 偏差改善更為明顯,24 h 預(yù)報(bào)時(shí)效內(nèi)偏差減小20%~50%,但同化地面PM10 觀測(cè)則是北方的預(yù)報(bào)效果提升更顯著;另外,對(duì)于京津冀地區(qū)而言,原模擬偏差隨預(yù)報(bào)時(shí)效在日變化上呈現(xiàn)正、負(fù)兩個(gè)方向的差異,顯示模式本身的排放、氣象和理化過程等多方面不確定性作用存在協(xié)同和制約作用,因而在某個(gè)時(shí)刻引入同化后初始場(chǎng)可能會(huì)導(dǎo)致后續(xù)某些時(shí)刻的偏差反而增大。

    (4)DA_1h 能夠在DA_6h 基礎(chǔ)上進(jìn)一步降低PM2.5 和PM10 的均方差并提高相關(guān)性,根據(jù)Sun et al.(2020)和本研究分別采用45 km 和9 km分辨率在中國地區(qū)的WRFDA-Chem 同化試驗(yàn),逐1 小時(shí)循環(huán)同化獲得的分析場(chǎng)能夠在不同程度上進(jìn)一步優(yōu)化預(yù)報(bào)效果并持續(xù)更長時(shí)間,對(duì)科學(xué)試驗(yàn)和構(gòu)建再分析數(shù)據(jù)集而言,采用逐1 小時(shí)同化是更好的選擇。但對(duì)于業(yè)務(wù)預(yù)報(bào)服務(wù)而言,需要同時(shí)考慮預(yù)報(bào)時(shí)效性和節(jié)約計(jì)算存儲(chǔ)資源,采用6 h 循環(huán)同化獲得的初始場(chǎng)啟動(dòng)預(yù)報(bào)已經(jīng)可以獲得較好的改善。

    由于云—?dú)馊苣z—輻射相互作用,提高PM2.5/PM10 預(yù)報(bào)準(zhǔn)確性能夠有效減小氣象要素的預(yù)報(bào)偏差,而對(duì)氣象要素預(yù)報(bào)的改進(jìn)也能反過來改進(jìn)PM2.5/PM10 預(yù)報(bào)。Yang et al.(2020)將RMAPSChem 實(shí)時(shí)模擬的氣溶膠濃度離線耦合到氣象預(yù)報(bào)模式RMAPS-ST,對(duì)中國華北地區(qū)進(jìn)行氣象預(yù)報(bào)試驗(yàn)發(fā)現(xiàn),采用實(shí)時(shí)模擬的氣溶膠信息比采用模式默認(rèn)AOD 固定值能夠更好地預(yù)報(bào)2 m 氣溫和10 m風(fēng)速,與不考慮氣溶膠反饋相比,實(shí)時(shí)更新氣溶膠信息能夠比采用AOD 固定值將2 m 氣溫和10 m風(fēng)速預(yù)報(bào)均方差分別進(jìn)一步較低2.8%和0.3%。王哲等(2014)利用WRF-NAQPMS 雙向耦合模式對(duì)京津冀秋季重霾的模擬表明,考慮氣溶膠—化學(xué)雙向反饋效應(yīng)能夠同時(shí)改善2 m 溫度、2 m 相對(duì)濕度、10 m 風(fēng)速和PM2.5 濃度模擬相關(guān)性和均方根誤差。本文聚焦同化在顆粒物污染預(yù)報(bào)中的作用,沒有考慮氣溶膠反饋的影響,但在后期研究中,將利用RMAPS-Chem 氣象—化學(xué)雙向反饋的模式框架與WRFDA-Chem 同化技術(shù)的應(yīng)用探討同化氣溶膠觀測(cè)資料對(duì)氣象預(yù)報(bào)和污染預(yù)報(bào)的影響。而在針對(duì)同化技術(shù)發(fā)展的未來工作中,一方面,需要對(duì)同化過程中觀測(cè)數(shù)據(jù)質(zhì)量控制環(huán)節(jié)進(jìn)行更多測(cè)試,以優(yōu)化同化效果;另一方面,同化預(yù)報(bào)效果的提升其實(shí)與模式本身性能息息相關(guān),需要在深入理解模式不確定性和偏差來源的情況下,進(jìn)一步開展模式和同化系統(tǒng)的協(xié)同發(fā)展。

    致謝感謝審稿專家和編輯為論文改進(jìn)提出的有益意見和建議,感謝中國環(huán)境監(jiān)測(cè)總站CNEMC 觀測(cè)網(wǎng)絡(luò)提供的數(shù)據(jù)。

    猜你喜歡
    氣溶膠偏差粒徑
    木屑粒徑對(duì)黑木耳栽培的影響試驗(yàn)*
    如何走出文章立意偏差的誤區(qū)
    氣溶膠傳播之謎
    兩矩形上的全偏差
    氣溶膠中210Po測(cè)定的不確定度評(píng)定
    基于近場(chǎng)散射的顆粒粒徑分布測(cè)量
    四川盆地秋季氣溶膠與云的相關(guān)分析
    關(guān)于均數(shù)與偏差
    Oslo結(jié)晶器晶體粒徑分布特征的CFD模擬
    大氣氣溶膠成核監(jiān)測(cè)
    亚洲欧美中文字幕日韩二区| 天天躁夜夜躁狠狠久久av| 亚洲国产av影院在线观看| 久久热在线av| 另类精品久久| 国产一区二区三区综合在线观看| 最黄视频免费看| 精品第一国产精品| 国产精品久久久久久精品古装| 久久这里有精品视频免费| 国产av精品麻豆| 日韩伦理黄色片| 国产一区二区 视频在线| 色婷婷久久久亚洲欧美| 国产又爽黄色视频| 精品国产一区二区三区久久久樱花| 久久精品国产自在天天线| 久久久久久久久免费视频了| 成人二区视频| 国产成人精品福利久久| 免费播放大片免费观看视频在线观看| 国产免费视频播放在线视频| 18禁观看日本| 久久久久国产网址| 亚洲美女视频黄频| 热re99久久国产66热| 飞空精品影院首页| 女人久久www免费人成看片| av免费观看日本| 久久久久久久久免费视频了| 18禁观看日本| 欧美最新免费一区二区三区| 亚洲精品乱久久久久久| 久久久久精品性色| 热99久久久久精品小说推荐| 国产精品成人在线| 美女xxoo啪啪120秒动态图| 97精品久久久久久久久久精品| 欧美97在线视频| 欧美xxⅹ黑人| 国产亚洲最大av| 狂野欧美激情性bbbbbb| 尾随美女入室| 国产精品无大码| 成人漫画全彩无遮挡| 一级a爱视频在线免费观看| 日韩视频在线欧美| 日韩av不卡免费在线播放| 午夜激情久久久久久久| 久久久久国产精品人妻一区二区| 亚洲欧美成人精品一区二区| 制服诱惑二区| 亚洲成人手机| 天天躁狠狠躁夜夜躁狠狠躁| 日韩人妻精品一区2区三区| 国产精品欧美亚洲77777| 美女国产高潮福利片在线看| kizo精华| 欧美bdsm另类| 嫩草影院入口| 亚洲国产看品久久| 久久99精品国语久久久| 一区二区三区四区激情视频| 新久久久久国产一级毛片| 免费高清在线观看日韩| 国产成人91sexporn| 久久午夜福利片| 97在线人人人人妻| 欧美精品一区二区大全| 18禁裸乳无遮挡动漫免费视频| 两个人看的免费小视频| 久久热在线av| 精品国产国语对白av| videossex国产| 久久久久精品性色| 欧美日韩亚洲高清精品| 久久久a久久爽久久v久久| 国产日韩欧美在线精品| 亚洲男人天堂网一区| 免费观看av网站的网址| 欧美日韩成人在线一区二区| 久久久久久伊人网av| 成人影院久久| 日韩免费高清中文字幕av| 三级国产精品片| 午夜福利视频在线观看免费| 免费在线观看完整版高清| 久久精品aⅴ一区二区三区四区 | 日本爱情动作片www.在线观看| 亚洲综合精品二区| 性色av一级| 精品亚洲成国产av| 欧美老熟妇乱子伦牲交| 欧美激情高清一区二区三区 | 日韩中字成人| 国产精品欧美亚洲77777| 王馨瑶露胸无遮挡在线观看| 精品少妇一区二区三区视频日本电影 | 男女免费视频国产| 欧美人与性动交α欧美精品济南到 | 亚洲欧洲日产国产| 嫩草影院入口| 国产无遮挡羞羞视频在线观看| 美女中出高潮动态图| 叶爱在线成人免费视频播放| av在线观看视频网站免费| 一本色道久久久久久精品综合| 综合色丁香网| 97人妻天天添夜夜摸| 亚洲国产精品国产精品| 欧美激情 高清一区二区三区| 日本-黄色视频高清免费观看| 国产熟女午夜一区二区三区| 国产亚洲午夜精品一区二区久久| 精品国产国语对白av| 妹子高潮喷水视频| 男女无遮挡免费网站观看| 日韩熟女老妇一区二区性免费视频| 亚洲精品视频女| 一级毛片 在线播放| 久久久国产一区二区| 亚洲,一卡二卡三卡| 午夜激情av网站| 老女人水多毛片| 国产欧美亚洲国产| 麻豆av在线久日| 精品久久蜜臀av无| 免费观看无遮挡的男女| 免费黄色在线免费观看| 精品亚洲成a人片在线观看| 欧美人与善性xxx| 成人毛片a级毛片在线播放| 一区二区av电影网| 久久热在线av| 欧美精品亚洲一区二区| 久久久国产欧美日韩av| 五月开心婷婷网| 亚洲人成电影观看| 一级片'在线观看视频| 热re99久久国产66热| 免费看av在线观看网站| 精品人妻在线不人妻| 叶爱在线成人免费视频播放| 欧美日韩精品成人综合77777| 男女边吃奶边做爰视频| av电影中文网址| 三上悠亚av全集在线观看| 国产不卡av网站在线观看| 叶爱在线成人免费视频播放| 又大又黄又爽视频免费| 欧美少妇被猛烈插入视频| 2021少妇久久久久久久久久久| 亚洲精品国产av蜜桃| 1024香蕉在线观看| av网站免费在线观看视频| av片东京热男人的天堂| 91在线精品国自产拍蜜月| 亚洲精品成人av观看孕妇| 国产又爽黄色视频| 精品国产露脸久久av麻豆| 中文字幕人妻熟女乱码| 伦理电影免费视频| 99久久精品国产国产毛片| 纯流量卡能插随身wifi吗| 人人妻人人爽人人添夜夜欢视频| a级毛片在线看网站| 亚洲美女视频黄频| 久久ye,这里只有精品| 亚洲av免费高清在线观看| 亚洲伊人色综图| 国产伦理片在线播放av一区| 欧美精品人与动牲交sv欧美| 日本vs欧美在线观看视频| 国产日韩一区二区三区精品不卡| 下体分泌物呈黄色| 欧美97在线视频| 国产有黄有色有爽视频| 亚洲av成人精品一二三区| 精品少妇久久久久久888优播| 一边摸一边做爽爽视频免费| 国产xxxxx性猛交| 狠狠婷婷综合久久久久久88av| 久久久久久伊人网av| 精品久久久精品久久久| 亚洲图色成人| 精品人妻熟女毛片av久久网站| 日产精品乱码卡一卡2卡三| 亚洲第一av免费看| 欧美老熟妇乱子伦牲交| 久久人妻熟女aⅴ| 免费在线观看黄色视频的| 日本av手机在线免费观看| 大片电影免费在线观看免费| 中文字幕最新亚洲高清| 久久97久久精品| 久久久精品区二区三区| 18禁国产床啪视频网站| 电影成人av| 日韩av免费高清视频| 中文字幕精品免费在线观看视频| 韩国精品一区二区三区| 国产精品一区二区在线观看99| 欧美日韩一级在线毛片| 人人澡人人妻人| 麻豆乱淫一区二区| 午夜免费男女啪啪视频观看| 亚洲精品自拍成人| 国产亚洲最大av| 国产乱来视频区| 黄片小视频在线播放| 99香蕉大伊视频| av女优亚洲男人天堂| 久久久久精品人妻al黑| 国产人伦9x9x在线观看 | 精品人妻偷拍中文字幕| 亚洲,一卡二卡三卡| 一区福利在线观看| 亚洲av成人精品一二三区| 亚洲欧美精品综合一区二区三区 | 久久影院123| 亚洲伊人色综图| 亚洲人成电影观看| 曰老女人黄片| 日韩一区二区视频免费看| 青春草国产在线视频| 人妻系列 视频| 精品99又大又爽又粗少妇毛片| 精品少妇黑人巨大在线播放| 波多野结衣一区麻豆| 校园人妻丝袜中文字幕| 中文字幕亚洲精品专区| 十八禁网站网址无遮挡| 成人黄色视频免费在线看| 在线观看美女被高潮喷水网站| 免费看av在线观看网站| 精品少妇黑人巨大在线播放| 欧美日韩视频精品一区| 精品少妇久久久久久888优播| 超碰97精品在线观看| 十八禁网站网址无遮挡| 男女边摸边吃奶| 国产精品三级大全| 777久久人妻少妇嫩草av网站| 日日啪夜夜爽| 你懂的网址亚洲精品在线观看| 一级毛片 在线播放| 久久影院123| 一级a爱视频在线免费观看| 午夜免费男女啪啪视频观看| 国产精品国产三级专区第一集| 亚洲欧洲精品一区二区精品久久久 | 三上悠亚av全集在线观看| 在线精品无人区一区二区三| 观看av在线不卡| 啦啦啦视频在线资源免费观看| 韩国精品一区二区三区| 国产精品.久久久| 99热国产这里只有精品6| 久久韩国三级中文字幕| 亚洲国产精品国产精品| 热99久久久久精品小说推荐| 国产黄色免费在线视频| 久久影院123| 国产精品久久久久久久久免| 久久久久久久久久久免费av| 日韩一卡2卡3卡4卡2021年| 午夜老司机福利剧场| 亚洲欧美色中文字幕在线| 91午夜精品亚洲一区二区三区| 男女边吃奶边做爰视频| 色播在线永久视频| 大码成人一级视频| 午夜免费鲁丝| 国产成人aa在线观看| 男人舔女人的私密视频| 成人午夜精彩视频在线观看| 黄频高清免费视频| 欧美国产精品一级二级三级| 熟女少妇亚洲综合色aaa.| 亚洲av男天堂| 啦啦啦在线免费观看视频4| 亚洲欧美日韩另类电影网站| 国产日韩欧美在线精品| 日韩人妻精品一区2区三区| 卡戴珊不雅视频在线播放| 中文字幕亚洲精品专区| 欧美+日韩+精品| 亚洲精品久久久久久婷婷小说| 少妇的丰满在线观看| 欧美精品av麻豆av| 国产日韩欧美视频二区| 十八禁网站网址无遮挡| 亚洲av.av天堂| 欧美精品亚洲一区二区| 在线观看美女被高潮喷水网站| 一区二区三区四区激情视频| 国产免费又黄又爽又色| 伊人亚洲综合成人网| 人人澡人人妻人| 999久久久国产精品视频| 亚洲欧美精品自产自拍| 美女xxoo啪啪120秒动态图| 亚洲精品中文字幕在线视频| 老司机影院毛片| 亚洲国产av影院在线观看| 久久久久国产精品人妻一区二区| 人人妻人人爽人人添夜夜欢视频| 最黄视频免费看| 久久国内精品自在自线图片| 精品一品国产午夜福利视频| 免费高清在线观看视频在线观看| 国产1区2区3区精品| 在线观看人妻少妇| 男女啪啪激烈高潮av片| 日本爱情动作片www.在线观看| 国产又爽黄色视频| 免费在线观看黄色视频的| 欧美 亚洲 国产 日韩一| 欧美激情 高清一区二区三区| 久久久久国产网址| 26uuu在线亚洲综合色| 欧美 亚洲 国产 日韩一| 亚洲欧美成人精品一区二区| 人人妻人人澡人人看| 五月天丁香电影| 考比视频在线观看| 日产精品乱码卡一卡2卡三| 国产成人免费无遮挡视频| 欧美激情 高清一区二区三区| 97在线视频观看| 免费少妇av软件| 国产精品亚洲av一区麻豆 | 波多野结衣av一区二区av| 国产麻豆69| 伊人久久大香线蕉亚洲五| 日本色播在线视频| 热re99久久精品国产66热6| 精品人妻偷拍中文字幕| 亚洲欧美日韩另类电影网站| 美女大奶头黄色视频| 亚洲av免费高清在线观看| 中文字幕另类日韩欧美亚洲嫩草| 一区二区三区乱码不卡18| 夫妻性生交免费视频一级片| 男女国产视频网站| 热99国产精品久久久久久7| 人妻 亚洲 视频| 性少妇av在线| 熟妇人妻不卡中文字幕| 国产极品粉嫩免费观看在线| 男人舔女人的私密视频| 日韩一卡2卡3卡4卡2021年| 亚洲图色成人| 成人国产av品久久久| 宅男免费午夜| 国产高清国产精品国产三级| 18禁国产床啪视频网站| 综合色丁香网| 美女高潮到喷水免费观看| 观看美女的网站| 91精品三级在线观看| 精品少妇久久久久久888优播| 亚洲成av片中文字幕在线观看 | 国产xxxxx性猛交| 七月丁香在线播放| 天堂俺去俺来也www色官网| 亚洲第一青青草原| 国产精品国产av在线观看| 在线观看免费日韩欧美大片| 嫩草影院入口| 午夜影院在线不卡| 日产精品乱码卡一卡2卡三| 免费观看无遮挡的男女| 成年女人在线观看亚洲视频| 久久精品国产自在天天线| a 毛片基地| 亚洲中文av在线| √禁漫天堂资源中文www| www.自偷自拍.com| 国产黄色免费在线视频| 美国免费a级毛片| 国产成人一区二区在线| 精品午夜福利在线看| 国产成人欧美| 国产女主播在线喷水免费视频网站| 一二三四中文在线观看免费高清| 男女无遮挡免费网站观看| 精品国产一区二区久久| 在线免费观看不下载黄p国产| 久久99一区二区三区| 国产精品熟女久久久久浪| 在线看a的网站| 亚洲欧美中文字幕日韩二区| 日产精品乱码卡一卡2卡三| 丰满饥渴人妻一区二区三| 午夜av观看不卡| 999精品在线视频| 欧美日韩精品成人综合77777| 精品99又大又爽又粗少妇毛片| 国产精品秋霞免费鲁丝片| 新久久久久国产一级毛片| 爱豆传媒免费全集在线观看| 午夜免费观看性视频| 成人影院久久| 免费大片黄手机在线观看| 亚洲成人手机| 国产免费福利视频在线观看| 最近2019中文字幕mv第一页| 两性夫妻黄色片| 女的被弄到高潮叫床怎么办| 日本vs欧美在线观看视频| 一个人免费看片子| 国产精品欧美亚洲77777| 18+在线观看网站| 青青草视频在线视频观看| 久久久久国产网址| 精品久久久精品久久久| 亚洲av欧美aⅴ国产| 又黄又粗又硬又大视频| 五月开心婷婷网| 欧美激情 高清一区二区三区| 日韩大片免费观看网站| av线在线观看网站| 男人操女人黄网站| 久久精品夜色国产| 中文字幕另类日韩欧美亚洲嫩草| 午夜福利视频精品| 老汉色av国产亚洲站长工具| 国产精品亚洲av一区麻豆 | 97在线视频观看| 十八禁高潮呻吟视频| 两性夫妻黄色片| 国产 精品1| 久久人妻熟女aⅴ| 亚洲成人手机| 国产在线视频一区二区| 亚洲av电影在线观看一区二区三区| 亚洲国产看品久久| 亚洲伊人久久精品综合| 丝袜在线中文字幕| 欧美日韩成人在线一区二区| 国产精品99久久99久久久不卡 | 天天操日日干夜夜撸| 国产在视频线精品| videosex国产| 欧美日韩视频高清一区二区三区二| 亚洲视频免费观看视频| 国产片内射在线| 亚洲成av片中文字幕在线观看 | 亚洲第一青青草原| av卡一久久| av.在线天堂| 久久久久国产精品人妻一区二区| 五月伊人婷婷丁香| 欧美精品高潮呻吟av久久| 一级片'在线观看视频| 宅男免费午夜| 午夜日韩欧美国产| av网站免费在线观看视频| 久久精品久久久久久久性| 日韩精品有码人妻一区| 欧美av亚洲av综合av国产av | 麻豆乱淫一区二区| 亚洲欧美成人综合另类久久久| 亚洲国产毛片av蜜桃av| 91久久精品国产一区二区三区| 精品少妇黑人巨大在线播放| 日韩免费高清中文字幕av| 天天躁夜夜躁狠狠久久av| 国产一级毛片在线| 十分钟在线观看高清视频www| 国产视频首页在线观看| h视频一区二区三区| 中文字幕最新亚洲高清| 性高湖久久久久久久久免费观看| 99久久精品国产国产毛片| 亚洲精品自拍成人| 男人舔女人的私密视频| videos熟女内射| 国产激情久久老熟女| 91精品三级在线观看| 午夜老司机福利剧场| a级毛片在线看网站| 如何舔出高潮| 久久久久视频综合| 亚洲第一青青草原| 咕卡用的链子| 丝袜人妻中文字幕| 日韩制服丝袜自拍偷拍| 国产精品三级大全| 黄色一级大片看看| 99久久精品国产国产毛片| 丝袜人妻中文字幕| 亚洲国产精品国产精品| 激情五月婷婷亚洲| 亚洲成色77777| 丝袜脚勾引网站| 亚洲美女黄色视频免费看| 精品人妻一区二区三区麻豆| 国产亚洲最大av| 国产白丝娇喘喷水9色精品| 成年人午夜在线观看视频| 久久99蜜桃精品久久| 老熟女久久久| av一本久久久久| av有码第一页| 久久久久网色| 中文字幕人妻熟女乱码| 国产色婷婷99| av国产精品久久久久影院| 狠狠婷婷综合久久久久久88av| 麻豆乱淫一区二区| 久久久久久伊人网av| 久久毛片免费看一区二区三区| 午夜福利视频精品| 男女啪啪激烈高潮av片| 黑人巨大精品欧美一区二区蜜桃| 亚洲精品国产一区二区精华液| 亚洲久久久国产精品| 精品久久久精品久久久| 大香蕉久久成人网| 新久久久久国产一级毛片| 国产人伦9x9x在线观看 | 丝袜脚勾引网站| 欧美日韩一级在线毛片| 成年动漫av网址| 日日摸夜夜添夜夜爱| 男女免费视频国产| 女人被躁到高潮嗷嗷叫费观| 2018国产大陆天天弄谢| 久久鲁丝午夜福利片| av不卡在线播放| 亚洲色图 男人天堂 中文字幕| 亚洲国产毛片av蜜桃av| 校园人妻丝袜中文字幕| 亚洲一码二码三码区别大吗| 中文字幕制服av| 有码 亚洲区| 日韩熟女老妇一区二区性免费视频| 人妻少妇偷人精品九色| 又大又黄又爽视频免费| 亚洲人成网站在线观看播放| 欧美日韩国产mv在线观看视频| 日韩一本色道免费dvd| av线在线观看网站| 中国国产av一级| 中文字幕另类日韩欧美亚洲嫩草| 亚洲成人一二三区av| 在线观看三级黄色| 久久久久精品人妻al黑| 亚洲欧美中文字幕日韩二区| 亚洲av成人精品一二三区| 国产午夜精品一二区理论片| 咕卡用的链子| 日本av免费视频播放| 日韩视频在线欧美| 91精品国产国语对白视频| 中文天堂在线官网| 亚洲国产看品久久| 久久韩国三级中文字幕| 国产黄色免费在线视频| 久久久国产欧美日韩av| 亚洲精品久久午夜乱码| 欧美精品人与动牲交sv欧美| 乱人伦中国视频| 久久影院123| 免费观看在线日韩| 女的被弄到高潮叫床怎么办| 国产精品 欧美亚洲| 亚洲经典国产精华液单| 青春草国产在线视频| 日韩中文字幕视频在线看片| 黄色配什么色好看| 侵犯人妻中文字幕一二三四区| 日本wwww免费看| 国产伦理片在线播放av一区| 老司机影院成人| 国产精品麻豆人妻色哟哟久久| 久久青草综合色| 边亲边吃奶的免费视频| 乱人伦中国视频| 国产精品女同一区二区软件| 亚洲一码二码三码区别大吗| www.自偷自拍.com| 久久热在线av| 永久免费av网站大全| 99久久人妻综合| 久久人人爽av亚洲精品天堂| 永久免费av网站大全| 亚洲精品一区蜜桃| 久久人人爽av亚洲精品天堂| 国产成人精品久久二区二区91 | 欧美97在线视频| av卡一久久| 久久精品亚洲av国产电影网| 国产成人精品一,二区| 亚洲精品av麻豆狂野| 婷婷色麻豆天堂久久| 久久久精品区二区三区| 大码成人一级视频| 久久久国产精品麻豆| 韩国高清视频一区二区三区| 午夜激情久久久久久久| 夫妻性生交免费视频一级片| 婷婷色av中文字幕| 欧美国产精品一级二级三级| 一级毛片黄色毛片免费观看视频| 女的被弄到高潮叫床怎么办| 国产亚洲最大av| 久久精品国产亚洲av高清一级| 丝瓜视频免费看黄片| 在线观看美女被高潮喷水网站| 美女大奶头黄色视频| 免费黄频网站在线观看国产| 日日摸夜夜添夜夜爱|