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

    基于INLA-SPDE方法的區(qū)域污染物模擬與預(yù)測(cè)

    2022-03-21 11:04:32澤,陳
    圖學(xué)學(xué)報(bào) 2022年1期
    關(guān)鍵詞:污染區(qū)域模型

    袁 澤,陳 斌

    基于INLA-SPDE方法的區(qū)域污染物模擬與預(yù)測(cè)

    袁 澤,陳 斌

    (北京大學(xué)地球與空間科學(xué)學(xué)院,北京 100871)

    采用傳統(tǒng)的空間插值方法對(duì)區(qū)域污染物進(jìn)行模擬與預(yù)測(cè),針對(duì)源數(shù)據(jù)分布不均,效果一般的問(wèn)題,提出了采用INLA-SPDE模型來(lái)模擬與預(yù)測(cè)區(qū)域污染物的方法。模型的空間分量使用隨機(jī)偏微分方程表達(dá),時(shí)間分量則采用一階時(shí)序自相關(guān)模型,同時(shí)還包含氣象參數(shù)等10種協(xié)變量,以2019年度京津冀地區(qū)日均PM2.5濃度為例,逐月建立了時(shí)空模擬與預(yù)測(cè)模型。實(shí)驗(yàn)結(jié)果表明,與經(jīng)典的克里金插值方法相比,在區(qū)域污染物分布的模擬上具有更好的效果,尤其在高值污染的預(yù)測(cè)上精度效果提升明顯,同時(shí)可得到區(qū)域污染風(fēng)險(xiǎn)等級(jí)等多種結(jié)果。進(jìn)一步基于模型的預(yù)測(cè)結(jié)果實(shí)現(xiàn)了京津冀地區(qū)日均PM2.5濃度時(shí)空可視化和虛擬仿真系統(tǒng),為普通民眾的出行或政府相關(guān)部門(mén)決策提供支持,驗(yàn)證了模型的實(shí)用性和價(jià)值。

    細(xì)顆粒物PM2.5;貝葉斯時(shí)空建模;INLA算法;仿真系統(tǒng);決策支持

    空氣污染一直是一個(gè)全球性問(wèn)題[1],隨著我國(guó)高人口密度的增加和快速工業(yè)化的結(jié)合,與經(jīng)濟(jì)增長(zhǎng)、工業(yè)化和城市化相關(guān)的嚴(yán)重空氣污染事件頻繁發(fā)生。當(dāng)下大氣顆粒物(particulate matter,PM),即懸浮在空氣當(dāng)中的固體顆?;蛞旱?,是最主要的污染物之一。而細(xì)顆粒物(空氣動(dòng)力學(xué)直徑小于等于2.5 μm的顆粒物稱(chēng)為PM2.5)因其粒徑更為微小,且具有最為明顯的污染效應(yīng)。同時(shí)許多流行病學(xué)和毒理學(xué)研究表明,PM2.5暴露會(huì)導(dǎo)致許多不良健康影響[2]。基于此,政府已經(jīng)將污染監(jiān)測(cè)和預(yù)警工作提到了一個(gè)新的高度。

    目前我國(guó)已經(jīng)建立了一定規(guī)模的地面空氣質(zhì)量監(jiān)測(cè)站點(diǎn)[3],其監(jiān)測(cè)結(jié)果可以有效地反映監(jiān)測(cè)站周邊的空氣質(zhì)量狀況。但受限于建設(shè)成本以及空間占地等問(wèn)題,地面監(jiān)測(cè)站數(shù)量還比較有限,且多位于發(fā)達(dá)大城市周邊,無(wú)法覆蓋大陸全域。同時(shí)近年來(lái)隨著監(jiān)測(cè)網(wǎng)絡(luò)的完善,產(chǎn)生了大量的時(shí)空數(shù)據(jù),這些數(shù)據(jù)大多具有較高的時(shí)間維度分辨率,一般為小時(shí)級(jí)甚至分鐘級(jí),而在空間維度上則相對(duì)稀疏,空間數(shù)據(jù)大多會(huì)隨著時(shí)間發(fā)生變化,對(duì)聯(lián)合時(shí)間與空間的數(shù)據(jù)模型的研究日益增多。

    在多種時(shí)空模型中,基于貝葉斯框架的時(shí)空建模方法近年來(lái)引起了越來(lái)越多學(xué)者的關(guān)注,其原因是時(shí)空數(shù)據(jù)具有內(nèi)在的不確定性,而貝葉斯方法從概率的角度出發(fā),能夠?qū)r(shí)空數(shù)據(jù)的不確定性進(jìn)行建模,最終給出概率上最優(yōu)的結(jié)果[4]。在研究空氣質(zhì)量時(shí)空分布模型時(shí),SARTO等[5]使用分層貝葉斯模型對(duì)意大利中部佩魯賈鎮(zhèn)的PM2.5濃度分布進(jìn)行了模擬,使用溫度、濕度等氣象因子作為協(xié)變量,同時(shí)添加空間效應(yīng)和時(shí)間自回歸項(xiàng),獲得了較好的擬合效果。

    目前對(duì)于貝葉斯時(shí)空模型參數(shù)的后驗(yàn)分布求解,大都采用馬爾可夫鏈蒙特卡洛方法(Markov chain Monto Carlo,MCMC),但由于到達(dá)馬爾可夫鏈的平穩(wěn)態(tài)往往需要很長(zhǎng)時(shí)間,在實(shí)際應(yīng)用中,MCMC方法并非作為一種常規(guī)分析工具。

    2009年,RUE等[6]針對(duì)一類(lèi)特殊的時(shí)空隱高斯隨機(jī)場(chǎng)的貝葉斯推斷,提出了一種稱(chēng)為集成嵌套拉普拉斯(integrated nested Laplace approximations,INLA)的快速求解算法。該方法可以求解一類(lèi)常見(jiàn)的結(jié)構(gòu)化加法回歸模型,這類(lèi)模型符合被一些超參數(shù)約束的隱高斯場(chǎng),同時(shí)具有非高斯響應(yīng)變量。通過(guò)使用集成的INLA逼近及其簡(jiǎn)化版本,該算法可以直接計(jì)算出后驗(yàn)分布非常精確的逼近,同時(shí)在MCMC算法需要運(yùn)行數(shù)小時(shí)或數(shù)天的情況下,該算法可在數(shù)秒或數(shù)分鐘內(nèi)提供更精確的估計(jì)。2011年,LINDGREN等[7]提出了一種隨機(jī)偏微分方程(stochastic partial differential equations,SPDE)算法,其可有效地鏈接連續(xù)高斯隨機(jī)場(chǎng)和離散的高斯馬爾科夫隨機(jī)場(chǎng)(Gaussian Markov random fields,GMRF),從而實(shí)現(xiàn)對(duì)基于Matern相關(guān)系數(shù)的連續(xù)分布的空間效應(yīng)的快速解算,以進(jìn)一步加速I(mǎi)NLA算法對(duì)貝葉斯時(shí)空模型后驗(yàn)分布的逼近。

    目前,該方法在國(guó)外已經(jīng)開(kāi)始普及,在與空氣質(zhì)量有關(guān)的研究方面,CAMELETTI等[8]針對(duì)意大利皮埃蒙特地區(qū)的健康狀況,采用INLA-SPDE方法獲得了區(qū)域污染物濃度,并開(kāi)展了對(duì)污染物濃度分布和區(qū)域健康風(fēng)險(xiǎn)關(guān)系的研究。FIORAVANTI等[9]基于意大利全國(guó)410個(gè)監(jiān)測(cè)站點(diǎn)2015年的日均PM10濃度數(shù)據(jù),使用INLA算法得到了意大利全域時(shí)空貝葉斯插值結(jié)果。但國(guó)內(nèi),該方法使用較少,幾乎只用于流行病時(shí)空分布的分析與模擬[10]。

    我國(guó)目前的空氣質(zhì)量地面監(jiān)測(cè)站主要分布在城市中心地帶,尚未形成空間分布均衡的監(jiān)測(cè)網(wǎng)絡(luò),因此難以對(duì)大陸區(qū)域空氣質(zhì)量狀況的時(shí)空分布進(jìn)行合理地評(píng)估。因此本文以京津冀地區(qū)為例,綜合利用研究區(qū)域內(nèi)的氣象數(shù)據(jù)和地理要素?cái)?shù)據(jù),建立了區(qū)域PM2.5濃度時(shí)空分布貝葉斯模型,并采取文獻(xiàn)[6-7]提出的INLA-SPDE加速算法,實(shí)現(xiàn)了對(duì)區(qū)域日均PM2.5濃度的快速模擬與預(yù)測(cè)。并以該預(yù)測(cè)結(jié)果為基礎(chǔ),構(gòu)建區(qū)域PM2.5濃度虛擬仿真決策支持系統(tǒng),實(shí)現(xiàn)京津冀區(qū)域PM2.5濃度時(shí)空分布圖、污染風(fēng)險(xiǎn)等級(jí)圖的可視化,支持與行政區(qū)劃圖、高程圖、人口圖等進(jìn)行疊加顯示,從而驗(yàn)證模型的結(jié)果實(shí)用性。該系統(tǒng)還支持交互放置特定的高污染源于地圖上,快速生成新情景下的空氣質(zhì)量空間分布圖,對(duì)比觀察污染物分布的變化,為污染物減排及控制等相關(guān)的空間問(wèn)題提供決策支持。

    1 原始數(shù)據(jù)處理與分析

    本文使用的原始PM2.5濃度數(shù)據(jù)來(lái)自京津冀地區(qū)236個(gè)國(guó)控空氣質(zhì)量監(jiān)測(cè)站,數(shù)據(jù)的原始時(shí)間分辨率為1 h,缺失值由站點(diǎn)當(dāng)前時(shí)刻前后3 h污染物濃度數(shù)據(jù)取均值代替。同時(shí)為了保證數(shù)據(jù)分布的正態(tài)性,去除了大于650的極端污染值數(shù)據(jù),日均濃度數(shù)據(jù)則由當(dāng)日24 h污染物濃度數(shù)據(jù)求平均值獲得,時(shí)間跨度從2019年1月到12月。最終得到月度PM2.5濃度日均值箱線圖(圖1),圖中的橘色虛線和紅色虛線分別對(duì)應(yīng)的PM2.5濃度為75 μg/m3和150 μg/m3,這是WHO標(biāo)準(zhǔn)中的第三級(jí)和第一級(jí)日均PM2.5濃度標(biāo)準(zhǔn),也是我國(guó)分級(jí)標(biāo)準(zhǔn)中的污染和重度污染分界線。

    圖1 京津冀地區(qū)2019年月度PM2.5濃度日均值箱線圖

    在2019年中,PM2.5濃度日均值在0~626 μg/m3之間,年均值為49.9 μg/m3,中位數(shù)為37 μg/m3,上下四分位數(shù)分別為22.2 μg/m3和60.8 μg/m3。該箱線圖與第2節(jié)研究一致,呈現(xiàn)出季節(jié)性規(guī)律:12月、1月、2月的PM2.5濃度日均值中位數(shù)大于50 μg/m3,夏季則小于30 μg/m3,其他月份在40 μg/m3。標(biāo)準(zhǔn)偏差也具有類(lèi)似的規(guī)律,1月和2月的標(biāo)準(zhǔn)偏差大于60 μg/m3,夏季小于20 μg/m3,其他月份在35 μg/m3左右。

    參與構(gòu)建貝葉斯模型的協(xié)變量在經(jīng)過(guò)反復(fù)試驗(yàn),最終選定了10個(gè)時(shí)空變量作為模型的預(yù)報(bào)因子。具體的數(shù)據(jù)信息見(jiàn)表1。為了避免數(shù)據(jù)自身的量級(jí)對(duì)研究的影響,所有的數(shù)據(jù)均進(jìn)行了標(biāo)準(zhǔn)化處理,即全部轉(zhuǎn)換為了均值為0,標(biāo)準(zhǔn)差為1的正態(tài)分布。

    表1 貝葉斯時(shí)空模型使用的協(xié)變量

    使用的6種氣象因子全部來(lái)自由美國(guó)國(guó)家環(huán)境中心(National Centers for Environment-Mental Prediction,NCEP)提供的6 h氣象再分析資料。該數(shù)據(jù)提供了每日0點(diǎn)、6點(diǎn)、12點(diǎn)、18點(diǎn)的空間分辨率為0.25度的格網(wǎng)化氣象要素,下載后以日為尺度重新計(jì)算了均值。氣溶膠光學(xué)厚度(aerosol optical depth,AOD)數(shù)據(jù)則來(lái)自歐空局(Copernicus Atmosphere Monitoring Service,CAMS)氣象再分析資料,是在空間分辨率為10 km的格網(wǎng)上對(duì)波長(zhǎng)為550 nm的AOD數(shù)值模擬得到的結(jié)果。

    京津冀地區(qū)各級(jí)道路數(shù)據(jù)獲取來(lái)自O(shè)SM (OpenStreetMap Project),從Geofabrik網(wǎng)站提供的服務(wù)下載得到。本文選取了OSM數(shù)據(jù)集中標(biāo)記為highway的線狀地物中的前四級(jí)標(biāo)簽作為基礎(chǔ)道路數(shù)據(jù),包括了京津冀地區(qū)的國(guó)道、省道以及縣道。最終投入模型的路網(wǎng)密度數(shù)據(jù)則是通過(guò)將京津冀地區(qū)劃分為10 km×10 km的格網(wǎng),計(jì)算每個(gè)格網(wǎng)內(nèi)道路數(shù)據(jù)的總長(zhǎng)度(km)來(lái)獲得的。海拔數(shù)據(jù)則來(lái)自USGS提供的GTOPO30全球數(shù)字高程模型,分辨率為30弧秒(約1 km)。人口數(shù)據(jù)采用WorldPop datasets提供的中國(guó)大陸1 km分辨率的人口密度數(shù)據(jù)。

    2 京津冀地區(qū)PM2.5日均濃度時(shí)空貝葉斯建模

    2.1 貝葉斯時(shí)空建模方案設(shè)計(jì)

    對(duì)2019年全年的京津冀地區(qū)日均PM2.5濃度進(jìn)行建模,按月份建立12個(gè)模型,具體模型如下:

    圖2 區(qū)域SPDE空間效應(yīng)三角網(wǎng)示意圖

    2.2 預(yù)測(cè)結(jié)果精度評(píng)估

    首先檢驗(yàn)?zāi)P偷恼军c(diǎn)預(yù)報(bào)精度,表2為12個(gè)月度模型在訓(xùn)練集和測(cè)試集上3次交叉驗(yàn)證后的平均結(jié)果,這些指標(biāo)數(shù)值均在轉(zhuǎn)換為原始數(shù)據(jù)尺度后計(jì)算獲得??傮w來(lái)說(shuō),12個(gè)模型在訓(xùn)練集和測(cè)試集上的預(yù)測(cè)均表現(xiàn)良好,而且效果相當(dāng),說(shuō)明模型并不存在過(guò)擬合。RMSE和MAE的值在冬季要高于夏季,這與冬季PM2.5濃度值整體偏高有關(guān)。

    表2 模型交叉驗(yàn)證結(jié)果(μg/m3)

    觀測(cè)值和預(yù)測(cè)值的關(guān)聯(lián)性如圖3所示,為了避免過(guò)多的散點(diǎn)圖,本文將12個(gè)月度模型按季度展示其真實(shí)值與預(yù)測(cè)值的一致性。圖中顏色越淺的部分代表處于此濃度的觀測(cè)數(shù)據(jù)的數(shù)量越多,紅色實(shí)線對(duì)應(yīng)1﹕1參考線。無(wú)論是在訓(xùn)練集還是測(cè)試集上,點(diǎn)幾乎都均勻地分布于參考線兩側(cè),同時(shí)淺色區(qū)域幾乎圍繞著參考線兩側(cè)最近處分布,可以認(rèn)為模型在4個(gè)季度對(duì)于大部分樣本數(shù)據(jù)預(yù)測(cè)效果表現(xiàn)良好。不過(guò)從圖3中也可以看到,在測(cè)試集上,當(dāng)觀測(cè)值大于250 μg/m3時(shí),出現(xiàn)了較多偏小的預(yù)測(cè)值。

    圖3 模型預(yù)測(cè)值與真實(shí)值對(duì)比圖((a)訓(xùn)練集;(b)測(cè)試集)

    圖4以污染最嚴(yán)重的1月和相對(duì)最輕的8月為例,繪制了區(qū)域內(nèi)236個(gè)站點(diǎn)月度的日均觀測(cè)值與預(yù)測(cè)置信區(qū)間的時(shí)序關(guān)聯(lián)圖,圖中黃色區(qū)域代表預(yù)測(cè)濃度值的95%置信區(qū)間。該時(shí)序圖表明京津冀地區(qū)PM2.5濃度日均值的時(shí)間變異性可以很好地被模型捕捉,觀測(cè)值濃度越高,預(yù)測(cè)值的置信區(qū)間范圍越大。但真實(shí)值總是落在預(yù)測(cè)的置信區(qū)間內(nèi)部,證明了模型對(duì)于時(shí)序特征的復(fù)現(xiàn)能力。

    2.3 預(yù)測(cè)結(jié)果對(duì)比分析與可視化

    本文選取具有代表性的日期,將模型的預(yù)測(cè)結(jié)果與傳統(tǒng)的空間插值通用方法——克里金插值方法進(jìn)行比較。即在相同的訓(xùn)練集上建立INLA-SPDE模型和進(jìn)行普通克里金空間插值,最終將測(cè)試集上的插值結(jié)果與INLA-SPDE模型在測(cè)試集上的預(yù)報(bào)結(jié)果進(jìn)行對(duì)比,其中貝葉斯模型的預(yù)測(cè)結(jié)果均采用模型的預(yù)測(cè)后驗(yàn)均值數(shù)據(jù)。

    圖4 PM2.5日均濃度時(shí)序驗(yàn)證圖((a)1月;(b)8月)

    具體選取的實(shí)驗(yàn)數(shù)據(jù)從2019年1月7日—18日,這一時(shí)段內(nèi)京津冀地區(qū)經(jīng)歷了一次污染從嚴(yán)重到消散的過(guò)程,具有一定的代表性。最終在完全相同的測(cè)試集上試驗(yàn)得到的RMSE數(shù)值見(jiàn)表3,可以發(fā)現(xiàn)使用的INLA-SPDE模型在空間上的預(yù)測(cè)效果要明顯優(yōu)于傳統(tǒng)的克里金插值。根據(jù)表3的數(shù)據(jù)繪制得到圖5,可以發(fā)現(xiàn)2種模型預(yù)測(cè)精度的差距在污染濃度越大時(shí)越明顯,證明了本文模型在高污染值上的預(yù)測(cè)能力遠(yuǎn)強(qiáng)于傳統(tǒng)的克里金插值方法。從PM2.5濃度預(yù)測(cè)的實(shí)際應(yīng)用出發(fā),對(duì)捕捉嚴(yán)重污染出現(xiàn)時(shí)機(jī)和位置的需求更大,因此本文模型具有較強(qiáng)的實(shí)用價(jià)值與現(xiàn)實(shí)意義。

    圖6以京津冀地區(qū)5 km×5 km的格網(wǎng)為基礎(chǔ),展示了克里金插值方法與本文模型的對(duì)比結(jié)果。選取了1月13日和15日分別為污染嚴(yán)重和開(kāi)始消散的代表日期,值得注意的是左右2幅圖顏色分級(jí)尺度不同。可以發(fā)現(xiàn),本文模型對(duì)于空間效應(yīng)的模擬與克里金插值相似,大尺度上的空間分布兩者保持著高度一致性。該模型對(duì)于污染濃度偏高的值更為敏感,能更為準(zhǔn)確地反映嚴(yán)重污染地帶的分布,2幅圖中極值亮點(diǎn)均出現(xiàn)在左圖,相比與傳統(tǒng)的克里金方法,INLA-SPDE模型對(duì)于小尺度上的空間變異性復(fù)現(xiàn)更為精確。

    表3 INLA模型與克里金插值預(yù)測(cè)RMSE結(jié)果(μg/m3)

    圖5 INLA模型與克里金插值預(yù)測(cè)RMSE對(duì)比圖

    圖6 INLA模型與克里金插值預(yù)測(cè)結(jié)果對(duì)比圖((a)1月13日;(b)1月15日)

    最后相較于傳統(tǒng)的克里金空間插值,因?yàn)樵撃P偷念A(yù)測(cè)結(jié)果本質(zhì)上是所有待預(yù)測(cè)格網(wǎng)點(diǎn)處污染濃度值的后驗(yàn)概率密度,因此還可以給出區(qū)域污染風(fēng)險(xiǎn)等級(jí)的空間分布,這對(duì)于人們?nèi)粘5某鲂邪才藕驼块T(mén)空間決策提供了支持。京津冀地區(qū)典型日期的污染風(fēng)險(xiǎn)如圖7所示。

    圖7 京津冀地區(qū)1月12日(左)和11月22日(右)污染風(fēng)險(xiǎn)等級(jí)空間分布圖

    3 區(qū)域PM2.5濃度虛擬仿真決策支持原型系統(tǒng)

    本節(jié)基于第2節(jié)建立的模型構(gòu)建京津冀地區(qū)PM2.5濃度虛擬仿真決策支持原型系統(tǒng)。主要目標(biāo)是驗(yàn)證該模型的實(shí)用性,實(shí)現(xiàn)任意時(shí)間模型結(jié)果的區(qū)域可視化,呈現(xiàn)時(shí)空預(yù)測(cè)結(jié)果中蘊(yùn)含的有效信息。并通過(guò)交互控制實(shí)現(xiàn)人工模擬大氣環(huán)境的改變來(lái)觀測(cè)區(qū)域效應(yīng),從而為普通民眾的出行或政府部門(mén)相關(guān)人員的決策提供支持。

    3.1 系統(tǒng)模塊設(shè)計(jì)

    原型系統(tǒng)包含3個(gè)基礎(chǔ)模塊,分別為模型層、數(shù)據(jù)層和交互層,如圖8所示。

    圖8 原型系統(tǒng)工作模塊圖

    將數(shù)據(jù)層中監(jiān)測(cè)站的實(shí)時(shí)數(shù)據(jù)與有關(guān)地理要素?cái)?shù)據(jù)合并成數(shù)據(jù)表,并存儲(chǔ)到Postgres數(shù)據(jù)庫(kù)中。氣象數(shù)據(jù)為格網(wǎng)化數(shù)據(jù),以GeoTIFF格式存儲(chǔ)。

    模型層輸入為選定日期的測(cè)站歷史數(shù)據(jù),調(diào)用相關(guān)月度的INLA-SPDE模型,將結(jié)果數(shù)據(jù)存儲(chǔ)到數(shù)據(jù)層后,再對(duì)區(qū)域PM2.5濃度預(yù)報(bào)的空間結(jié)果可視化。

    交互模塊通過(guò)Web前端界面加入相應(yīng)交互按鈕,支持人工模擬新污染源的添加與修改,求解得到交互后的模擬空氣質(zhì)量時(shí)空分布結(jié)果并可視化。

    3.2 模型的原型系統(tǒng)驗(yàn)證

    用于驗(yàn)證的原型系統(tǒng)基本界面主要包括,交互和結(jié)果展示2個(gè)區(qū)域。交互區(qū)如圖9左側(cè)邊欄所示,包含日期選擇、想要修改污染值的點(diǎn)位經(jīng)緯度信息、污染輻射范圍和濃度值的錄入,其中位置的經(jīng)緯度信息也可以通過(guò)地圖區(qū)的繪制功能交互點(diǎn)選實(shí)現(xiàn)。

    圖9為主區(qū)域?yàn)榫┙蚣降貐^(qū)PM2.5濃度時(shí)空分布結(jié)果地圖,系統(tǒng)可根據(jù)所選可視化區(qū)域PM2.5濃度空間分布模擬結(jié)果。地圖區(qū)通過(guò)左側(cè)的功能欄,支持顯示污染風(fēng)險(xiǎn)等級(jí)圖和空間分布結(jié)果圖與多種驗(yàn)證數(shù)據(jù)圖的疊加。該結(jié)果展示區(qū)還包括若干功能,如上方顯示的是PM2.5濃度超標(biāo)的測(cè)站個(gè)數(shù)、預(yù)測(cè)PM2.5濃度超標(biāo)的區(qū)域占京津冀地區(qū)的比例和重污染風(fēng)險(xiǎn)濃度在75%以上的區(qū)域比例,對(duì)于整個(gè)區(qū)域這一時(shí)刻的污染情況可以有直觀且清晰的感知。右側(cè)則顯示了污染具體的行政區(qū)劃分布情況,包括站點(diǎn)PM2.5濃度排序和以市為基本單位計(jì)算的高污染風(fēng)險(xiǎn)區(qū)域所占的比例,從而為特定城市的相關(guān)人群提供有價(jià)值的參考信息。

    界面下方為交互模擬結(jié)果展示,根據(jù)交互區(qū)輸入的相關(guān)信息得到新的空間分布結(jié)果并可視化,與主地圖區(qū)對(duì)比,效果如圖10所示。圖中更改濃度值的點(diǎn)位為原結(jié)果圖中的藍(lán)色圓圈,將此處的污染值向下修正到了50 μg/m3,并設(shè)定輻射范圍為100 km。在變動(dòng)結(jié)果圖中,可以清楚地看到藍(lán)色圓圈周?chē)谋O(jiān)測(cè)站值全部降低到75 μg/m3以下(圖中的紅色代表污染值在75 μg/m3以上的站點(diǎn),綠色代表在75 μg/m3以下的站點(diǎn))。從區(qū)域分布上可知,整個(gè)北京以北的地區(qū)污染狀況顯著改善,為了量化這種污染改善的程度,系統(tǒng)右下角還展示了不同污染等級(jí)的區(qū)域格網(wǎng)數(shù)目變動(dòng)百分比,可以看到重度污染以上的區(qū)域數(shù)目大幅減少,而轉(zhuǎn)為良好等級(jí)的區(qū)域數(shù)目則明顯上升。

    圖9 原型系統(tǒng)基礎(chǔ)界面圖

    (a) (b)

    因此,通過(guò)該交互控制功能,可以清楚地展現(xiàn)在特定區(qū)域減排或設(shè)立新的污染源的影響程度,為相關(guān)政策的針對(duì)性制定和部分建筑的空間選址問(wèn)題提供了重要的參考,真正實(shí)現(xiàn)了空間決策與支持功能,進(jìn)一步驗(yàn)證了模型的實(shí)用價(jià)值。

    4 結(jié)論與反思

    本文提出了基于INLA-SPDE模型的區(qū)域污染物模擬與預(yù)測(cè)方法,實(shí)現(xiàn)了京津冀地區(qū)PM2.5濃度日均值的時(shí)空建模與預(yù)測(cè),驗(yàn)證了模型預(yù)測(cè)精度,得到了京津冀地區(qū)5 km×5 km格網(wǎng)的日均PM2.5濃度時(shí)空分布預(yù)報(bào)圖和污染風(fēng)險(xiǎn)等級(jí)圖等。結(jié)果表明,INLA-SPDE方法在區(qū)域污染模擬與預(yù)測(cè)上相比傳統(tǒng)空間插值方法具有更好的效果,尤其對(duì)于高污染值的預(yù)測(cè)更為精確,并能給出區(qū)域污染風(fēng)險(xiǎn)等級(jí)預(yù)測(cè)結(jié)果,在區(qū)域污染的預(yù)警和防控上具有重要的現(xiàn)實(shí)意義。基于該模型構(gòu)建了京津冀地區(qū)PM2.5濃度預(yù)報(bào)與控制決策支持原型系統(tǒng),解決了區(qū)域級(jí)別的PM2.5濃度時(shí)空可視化的問(wèn)題。通過(guò)將人工模擬得到的仿真結(jié)果與原始結(jié)果對(duì)比,為人類(lèi)的空間決策行為提供了有效地支持,驗(yàn)證了該模型的實(shí)用性與價(jià)值。

    [1] Report from Intergovernmental Panel on Climate Change (IPCC). Climate change 2007: the physical science basis[EB/OL]. [2021-04-01]. http://www.doc88.com/p- 293361412559.html.

    [2] APTE J S, BRAUER M, COHEN A J, et al. Ambient PM2.5reduces global and regional life expectancy[J]. Environmental Science & Technology Letters, 2018, 5(9): 546-551.

    [3] CAI S Y, MA Q, WANG S X, et al. Impact of air pollution control policies on future PM2.5concentrations and their source contributions in China[J]. Journal of Environmental Management, 2018, 227: 124-133.

    [4] GAUDARD M, KARSON M, LINDER E, et al. Bayesian spatial prediction[J]. Environmental and Ecological Statistics, 1999, 6(2): 147-171.

    [5] SARTO S D, RANALLI M G, BAKAR K S, et al. Bayesian spatiotemporal modeling of urban air pollution dynamics[M]// Topics on Methodological and Applied Statistical Inference. Cham: Springer International Publishing, 2016: 95-103.

    [6] RUE H, MARTINO S, CHOPIN N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[J]. Journal of the Royal Statistical Society: Series B: Statistical Methodology, 2009, 71(2): 319-392.

    [7] LINDGREN F, RUE H, LINDSTR?M J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[J]. Journal of the Royal Statistical Society: Series B: Statistical Methodology, 2011, 73(4): 423-498.

    [8] CAMELETTI M, GóMEZ-RUBIO V, BLANGIARDO M. Bayesian modelling for spatially misaligned health and air pollution data through the INLA-SPDE approach[J]. Spatial Statistics, 2019, 31: 100353.

    [9] FIORAVANTI G, MARTINO S, CAMELETTI M, et al. Spatio-temporal modelling of PM10daily concentrations in Italy using the SPDE approach[J]. Atmospheric Environment, 2021, 248: 118192.

    [10] 丁書(shū)珍, 張輝國(guó), 胡錫健, 等. 利用R-INLA方法分析宏觀因素對(duì)艾滋病疫情的影響[J]. 中國(guó)艾滋病性病, 2018, 24(12): 1192-1196.

    DING S Z, ZHANG H G, HU X J, et al. Influence of macroscopic factors on HIV/AIDS epidemic based on R-INLA method[J]. Chinese Journal of AIDS & STD, 2018, 24(12): 1192-1196 (in Chinese).

    [11] BLANGIARDO M, CAMELETTI M. Spatial and spatio-temporal Bayesian models with R-INLA[M]. Chichester: John Wiley & Sons, Ltd, 2015: 268-283.

    [12] SIMPSON D, RUE H, RIEBLER A, et al. Penalising model component complexity: a principled, practical approach to constructing priors[EB/OL]. [2021-04-04]. https://arxiv.org/ pdf/1403.4630.pdf.

    [13] S?RBYE S H, RUE H. Penalised complexity priors for stationary autoregressive processes[J]. Journal of Time Series Analysis, 2017, 38(6): 923-935.

    Simulation and prediction of regional pollutants based on INLA-SPDE method

    YUAN Ze, CHEN Bin

    (School of Earth and Space Sciences, Peking University, Beijing 100871, China)

    The simulation and prediction of regional pollutants generally use the traditional spatial interpolation method, which cannot obtain accurate results when the source data is not uniformly distributed. To address these problems, a method for simulation and prediction of regional pollutants based on the INLA-SPDE model was proposed. The interpolation model was based on a Bayesian hierarchical model where the spatial-component was represented through the stochastic partial differential equation (SPDE) approach, with a lag-1 temporal autoregressive component (AR1). In addition, the model included 10 spatial and spatio-temporal predictors such as meteorological variables. By building 12 models for each month with the integrated nested Laplace approximation (INLA), this research realized the spatio-temporal simulation and prediction of PM2.5concentration at daily resolution in the Beijing-Tianjin-Hebei region in 2019. Experiments show that compared with traditional Kriging interpolation methods, the proposed model can yield a better prediction of air pollutants at regional scale. Particularly, the prediction accuracy of high-value pollutants was improved significantly, and air pollutants exceedance probabilities can also be generated. Furthermore, a system for regional PM2.5concentration simulation and decision support was established, the system can provide support for the travel of ordinary people or the decision-making of government officials, and verify the practicability and value of the proposed model.

    PM2.5; spatio-temporal Bayesian hierarchical model; integrated nested Laplace approximation; simulation system; decision support

    18 May,2021;

    TP 391

    10.11996/JG.j.2095-302X.2022010125

    A

    2095-302X(2022)01-0125-08

    2021-05-18;

    2021-08-21

    21 August,2021

    袁 澤(1992–),男,碩士研究生。主要研究方向?yàn)樘摂M地理環(huán)境。E-mail:180121089@pku.edu.cn

    YUAN Ze (1992–), master student. His main research interest covers virtual geographic environments. E-mail:180121089@pku.edu.cn

    陳 斌(1973–),男,教授,博士。主要研究方向?yàn)樘摂M地理環(huán)境。E-mail:gischen@pku.edu.cn

    CHEN Bin (1973–), professor, Ph.D. His main research interest covers virtual geographic environments.E-mail:gischen@pku.edu.cn

    猜你喜歡
    污染區(qū)域模型
    一半模型
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    堅(jiān)決打好污染防治攻堅(jiān)戰(zhàn)
    堅(jiān)決打好污染防治攻堅(jiān)戰(zhàn)
    3D打印中的模型分割與打包
    關(guān)于四色猜想
    分區(qū)域
    基于嚴(yán)重區(qū)域的多PCC點(diǎn)暫降頻次估計(jì)
    對(duì)抗塵污染,遠(yuǎn)離“霾”伏
    都市麗人(2015年5期)2015-03-20 13:33:49
    成人性生交大片免费视频hd| 亚州av有码| 色播亚洲综合网| 99热6这里只有精品| 国产色婷婷99| 日本-黄色视频高清免费观看| 99久久精品热视频| 女人十人毛片免费观看3o分钟| 在线免费观看不下载黄p国产| 天堂网av新在线| 亚洲av二区三区四区| 美女大奶头视频| 国产一区二区亚洲精品在线观看| 久久久久精品性色| 日韩视频在线欧美| 国产午夜福利久久久久久| 菩萨蛮人人尽说江南好唐韦庄| 国产黄片视频在线免费观看| 嘟嘟电影网在线观看| 亚洲人成网站在线观看播放| av在线蜜桃| 亚洲综合色惰| 国模一区二区三区四区视频| 久久久a久久爽久久v久久| 欧美最新免费一区二区三区| 亚洲av男天堂| 午夜老司机福利剧场| 高清毛片免费看| 九草在线视频观看| 美女大奶头视频| 欧美日韩视频高清一区二区三区二| 少妇人妻精品综合一区二区| 午夜免费男女啪啪视频观看| freevideosex欧美| 18禁在线播放成人免费| 免费看日本二区| 亚洲,欧美,日韩| 日本与韩国留学比较| 激情 狠狠 欧美| 少妇人妻精品综合一区二区| 美女主播在线视频| 一级二级三级毛片免费看| 久久精品夜色国产| 在线播放无遮挡| 婷婷色麻豆天堂久久| 97在线视频观看| 亚洲欧美中文字幕日韩二区| 亚洲国产日韩欧美精品在线观看| 免费看a级黄色片| 成人漫画全彩无遮挡| 国产精品国产三级国产专区5o| 午夜视频国产福利| 免费高清在线观看视频在线观看| 久久精品久久精品一区二区三区| 婷婷色av中文字幕| 精品酒店卫生间| 欧美bdsm另类| 精品久久久久久久末码| 又黄又爽又刺激的免费视频.| 精品国产一区二区三区久久久樱花 | 中文字幕亚洲精品专区| 看十八女毛片水多多多| 精品久久久久久电影网| 毛片女人毛片| 久久这里有精品视频免费| 国产永久视频网站| 欧美性猛交╳xxx乱大交人| 麻豆乱淫一区二区| 人妻一区二区av| 国产午夜精品久久久久久一区二区三区| 青春草亚洲视频在线观看| 少妇熟女欧美另类| 人妻制服诱惑在线中文字幕| 亚洲av男天堂| 好男人在线观看高清免费视频| 夜夜爽夜夜爽视频| 国内少妇人妻偷人精品xxx网站| 国产黄a三级三级三级人| 亚洲国产色片| 国产成人一区二区在线| 26uuu在线亚洲综合色| 亚洲aⅴ乱码一区二区在线播放| 亚洲精品国产av蜜桃| 国产成人精品福利久久| 国产精品久久久久久精品电影| 午夜福利在线观看吧| 日韩欧美国产在线观看| 国产精品一区二区性色av| 一边亲一边摸免费视频| 日韩 亚洲 欧美在线| 赤兔流量卡办理| 亚洲精品,欧美精品| 嫩草影院新地址| 女的被弄到高潮叫床怎么办| 一级毛片黄色毛片免费观看视频| 国产精品日韩av在线免费观看| 一级黄片播放器| 国产成人一区二区在线| 最新中文字幕久久久久| 男插女下体视频免费在线播放| 亚洲人成网站高清观看| 久久精品人妻少妇| 日本av手机在线免费观看| 精品99又大又爽又粗少妇毛片| 成年版毛片免费区| 国产成年人精品一区二区| 亚洲国产色片| 搞女人的毛片| 亚洲精品乱久久久久久| 国产中年淑女户外野战色| 亚洲精品一二三| 国产一区二区三区av在线| 成人午夜高清在线视频| 久久久久久久亚洲中文字幕| 日产精品乱码卡一卡2卡三| 日本av手机在线免费观看| 中文字幕免费在线视频6| 一个人看的www免费观看视频| 我的女老师完整版在线观看| 久久久久精品性色| 国内精品一区二区在线观看| 久久久国产一区二区| 在线免费观看不下载黄p国产| 欧美+日韩+精品| 又大又黄又爽视频免费| 一本一本综合久久| 少妇的逼水好多| 人人妻人人澡欧美一区二区| 午夜免费激情av| 亚洲丝袜综合中文字幕| 精品久久久久久久久亚洲| 能在线免费看毛片的网站| 欧美97在线视频| 婷婷色麻豆天堂久久| 国产高清有码在线观看视频| 久久久久久久亚洲中文字幕| 日韩一本色道免费dvd| 欧美性感艳星| 高清毛片免费看| 大话2 男鬼变身卡| 一级二级三级毛片免费看| 亚洲久久久久久中文字幕| 肉色欧美久久久久久久蜜桃 | 国产精品一区www在线观看| 欧美高清成人免费视频www| 97人妻精品一区二区三区麻豆| 亚洲高清免费不卡视频| 一级毛片aaaaaa免费看小| 自拍偷自拍亚洲精品老妇| 卡戴珊不雅视频在线播放| 国产色爽女视频免费观看| 亚洲精品乱码久久久久久按摩| 欧美日韩视频高清一区二区三区二| 亚洲图色成人| 纵有疾风起免费观看全集完整版 | 成人综合一区亚洲| 国产真实伦视频高清在线观看| 亚洲婷婷狠狠爱综合网| 秋霞伦理黄片| 少妇被粗大猛烈的视频| 成人av在线播放网站| 亚洲性久久影院| 亚洲美女视频黄频| 51国产日韩欧美| 国产一区二区三区综合在线观看 | 免费播放大片免费观看视频在线观看| 麻豆成人午夜福利视频| 国产一区二区三区综合在线观看 | 亚洲av成人精品一二三区| 99久国产av精品国产电影| kizo精华| 色综合亚洲欧美另类图片| 简卡轻食公司| 久久久久国产网址| 天堂√8在线中文| 啦啦啦啦在线视频资源| 免费观看av网站的网址| 日韩成人av中文字幕在线观看| av线在线观看网站| 久久久亚洲精品成人影院| 一夜夜www| 十八禁国产超污无遮挡网站| 日日啪夜夜撸| 美女cb高潮喷水在线观看| 久久久久久久久中文| eeuss影院久久| 日韩亚洲欧美综合| 天天一区二区日本电影三级| 国产成人福利小说| 亚洲久久久久久中文字幕| 国产女主播在线喷水免费视频网站 | 国产一区有黄有色的免费视频 | 国产老妇伦熟女老妇高清| 国产亚洲精品av在线| 亚洲国产av新网站| 久久久久性生活片| 精品熟女少妇av免费看| 久久99热这里只频精品6学生| 国产综合懂色| 91av网一区二区| 亚洲精品国产av成人精品| 国产在线男女| 久久久久久久久久黄片| 国产精品伦人一区二区| 久久久久久九九精品二区国产| 啦啦啦中文免费视频观看日本| 亚洲人成网站高清观看| 国产一区二区在线观看日韩| 三级毛片av免费| 成年版毛片免费区| 最近视频中文字幕2019在线8| 丝瓜视频免费看黄片| 国内少妇人妻偷人精品xxx网站| 亚洲aⅴ乱码一区二区在线播放| 国产探花在线观看一区二区| 99久久精品热视频| 亚洲成人一二三区av| 在线天堂最新版资源| 搡老乐熟女国产| 天堂影院成人在线观看| 91av网一区二区| 成年人午夜在线观看视频 | 国产精品爽爽va在线观看网站| 亚洲美女搞黄在线观看| 只有这里有精品99| 日韩亚洲欧美综合| 精品国产一区二区三区久久久樱花 | 国产精品一区二区三区四区免费观看| 天堂√8在线中文| 成人亚洲精品av一区二区| av网站免费在线观看视频 | 欧美+日韩+精品| 亚洲精品色激情综合| 国产单亲对白刺激| 日韩av不卡免费在线播放| 欧美成人一区二区免费高清观看| 亚洲在久久综合| 日韩大片免费观看网站| 嘟嘟电影网在线观看| 一区二区三区高清视频在线| 97热精品久久久久久| 久久久久免费精品人妻一区二区| 2018国产大陆天天弄谢| 亚洲精品乱码久久久v下载方式| 男女边摸边吃奶| 国模一区二区三区四区视频| 十八禁网站网址无遮挡 | 特级一级黄色大片| 国产精品日韩av在线免费观看| 日日啪夜夜爽| 少妇丰满av| 亚洲欧美一区二区三区国产| 三级国产精品欧美在线观看| 能在线免费看毛片的网站| 欧美+日韩+精品| 九色成人免费人妻av| 日本三级黄在线观看| 国产精品国产三级专区第一集| 亚洲欧美清纯卡通| 人人妻人人澡人人爽人人夜夜 | 一级毛片aaaaaa免费看小| 日本wwww免费看| 久久人人爽人人片av| 亚洲第一区二区三区不卡| 69人妻影院| 久久这里有精品视频免费| 亚洲av.av天堂| 嘟嘟电影网在线观看| 日韩一本色道免费dvd| 久久久久九九精品影院| 久久久久久久久久人人人人人人| 成人午夜高清在线视频| 日韩国内少妇激情av| 国产黄a三级三级三级人| 日日啪夜夜爽| 免费观看a级毛片全部| 国产麻豆成人av免费视频| 插逼视频在线观看| 国产精品久久视频播放| 久久久久久久久久久免费av| 国产黄片视频在线免费观看| 91狼人影院| 免费不卡的大黄色大毛片视频在线观看 | 精品一区二区三区人妻视频| 精品久久久精品久久久| 国产男女超爽视频在线观看| 大片免费播放器 马上看| 亚洲性久久影院| 夫妻午夜视频| 美女被艹到高潮喷水动态| 日日撸夜夜添| 天堂影院成人在线观看| 久久久成人免费电影| 亚洲一区高清亚洲精品| 免费高清在线观看视频在线观看| 久久久久精品性色| 国产免费一级a男人的天堂| 美女内射精品一级片tv| 国产毛片a区久久久久| 亚洲av电影在线观看一区二区三区 | 97超碰精品成人国产| 在线免费观看的www视频| 18禁动态无遮挡网站| 婷婷色麻豆天堂久久| 亚洲国产精品成人久久小说| 大话2 男鬼变身卡| 亚洲国产最新在线播放| 老司机影院毛片| 久久97久久精品| 国产午夜精品论理片| 夫妻午夜视频| 精品久久久久久久人妻蜜臀av| 亚洲av免费在线观看| 亚洲美女搞黄在线观看| 亚洲av成人精品一二三区| 国产精品久久久久久久久免| 97人妻精品一区二区三区麻豆| 日韩一区二区三区影片| 在线观看av片永久免费下载| 在线免费十八禁| 午夜福利在线观看免费完整高清在| 大片免费播放器 马上看| 人妻制服诱惑在线中文字幕| 国产av国产精品国产| 久久久久免费精品人妻一区二区| 久久综合国产亚洲精品| 亚洲,欧美,日韩| 最近视频中文字幕2019在线8| 国产伦精品一区二区三区四那| 777米奇影视久久| 久久久a久久爽久久v久久| 亚洲伊人久久精品综合| 亚洲天堂国产精品一区在线| 国产一区二区三区av在线| 日本欧美国产在线视频| 亚洲欧美成人精品一区二区| 免费观看av网站的网址| 一级av片app| 十八禁网站网址无遮挡 | 欧美性猛交╳xxx乱大交人| 建设人人有责人人尽责人人享有的 | 久久久久久伊人网av| 51国产日韩欧美| 三级男女做爰猛烈吃奶摸视频| 欧美性感艳星| 久久精品熟女亚洲av麻豆精品 | 亚洲人与动物交配视频| 国产精品久久久久久久电影| 嫩草影院新地址| 波多野结衣巨乳人妻| 丰满人妻一区二区三区视频av| 91精品国产九色| 日韩欧美 国产精品| av在线播放精品| 免费观看精品视频网站| 黄片wwwwww| 中文乱码字字幕精品一区二区三区 | 国产亚洲午夜精品一区二区久久 | 久久久亚洲精品成人影院| 国产人妻一区二区三区在| 两个人的视频大全免费| 精品人妻熟女av久视频| 日韩av不卡免费在线播放| 亚洲欧美一区二区三区国产| 蜜臀久久99精品久久宅男| 麻豆乱淫一区二区| 欧美日韩在线观看h| 干丝袜人妻中文字幕| 亚洲欧美精品专区久久| 免费看日本二区| 亚洲熟女精品中文字幕| 亚洲精品中文字幕在线视频 | 人妻系列 视频| 亚洲av中文av极速乱| 97精品久久久久久久久久精品| 日韩av不卡免费在线播放| 狂野欧美白嫩少妇大欣赏| 男人舔奶头视频| 亚洲精品第二区| 91午夜精品亚洲一区二区三区| 国产美女午夜福利| 亚洲精品久久久久久婷婷小说| 三级毛片av免费| 国产黄色免费在线视频| 久久精品久久久久久久性| 街头女战士在线观看网站| 久久99热6这里只有精品| 91aial.com中文字幕在线观看| 免费看日本二区| 日韩中字成人| 又爽又黄a免费视频| 七月丁香在线播放| 欧美一级a爱片免费观看看| 中文精品一卡2卡3卡4更新| 夜夜爽夜夜爽视频| 91久久精品国产一区二区三区| 国产精品不卡视频一区二区| 日本爱情动作片www.在线观看| 看免费成人av毛片| 亚洲婷婷狠狠爱综合网| 寂寞人妻少妇视频99o| av女优亚洲男人天堂| 22中文网久久字幕| 国产亚洲一区二区精品| 看非洲黑人一级黄片| 午夜福利成人在线免费观看| 亚洲伊人久久精品综合| 黄色欧美视频在线观看| 狂野欧美激情性xxxx在线观看| av在线亚洲专区| 啦啦啦中文免费视频观看日本| 成年版毛片免费区| 99久久精品国产国产毛片| 精品一区二区三区人妻视频| 99热全是精品| 日韩伦理黄色片| 亚洲美女搞黄在线观看| 国产成年人精品一区二区| 国产淫语在线视频| 我要看日韩黄色一级片| 亚洲欧美日韩东京热| 大片免费播放器 马上看| 一级毛片黄色毛片免费观看视频| 国产精品一区www在线观看| 2022亚洲国产成人精品| 国产v大片淫在线免费观看| 搡女人真爽免费视频火全软件| 中文字幕人妻熟人妻熟丝袜美| 国产精品人妻久久久影院| 亚洲av中文字字幕乱码综合| 日本av手机在线免费观看| 成年版毛片免费区| 2021少妇久久久久久久久久久| 久久久久久九九精品二区国产| 三级国产精品片| 深夜a级毛片| 中文字幕免费在线视频6| 国产成人免费观看mmmm| 可以在线观看毛片的网站| or卡值多少钱| 伦理电影大哥的女人| 精品久久久精品久久久| 精品人妻偷拍中文字幕| 女人被狂操c到高潮| 91av网一区二区| 色尼玛亚洲综合影院| 爱豆传媒免费全集在线观看| 国产成人午夜福利电影在线观看| 亚洲欧美精品自产自拍| 精品一区二区三区视频在线| 国产精品福利在线免费观看| 九九爱精品视频在线观看| 最近最新中文字幕免费大全7| 啦啦啦韩国在线观看视频| 舔av片在线| 看免费成人av毛片| 99久久精品国产国产毛片| 男女边摸边吃奶| 国产v大片淫在线免费观看| 久久久久久伊人网av| 日韩电影二区| 搡老妇女老女人老熟妇| 国产有黄有色有爽视频| 乱系列少妇在线播放| 人人妻人人澡欧美一区二区| 99久久人妻综合| 国产av不卡久久| 日韩在线高清观看一区二区三区| 一二三四中文在线观看免费高清| 成年女人看的毛片在线观看| 秋霞在线观看毛片| 亚洲国产精品sss在线观看| 欧美成人午夜免费资源| 国产成人精品婷婷| 国产伦理片在线播放av一区| 欧美xxxx黑人xx丫x性爽| 中文天堂在线官网| 蜜桃亚洲精品一区二区三区| 国产 亚洲一区二区三区 | 免费少妇av软件| 国产精品爽爽va在线观看网站| 欧美激情久久久久久爽电影| 2021少妇久久久久久久久久久| freevideosex欧美| 99久久人妻综合| 日本猛色少妇xxxxx猛交久久| 久久精品国产亚洲av涩爱| 精品久久久精品久久久| 神马国产精品三级电影在线观看| 亚洲久久久久久中文字幕| 国产亚洲91精品色在线| 久久精品国产亚洲av涩爱| 一级毛片我不卡| 亚洲激情五月婷婷啪啪| 免费黄网站久久成人精品| 国产精品一区二区三区四区久久| 91狼人影院| 成人无遮挡网站| 天天躁夜夜躁狠狠久久av| 国产精品蜜桃在线观看| 免费不卡的大黄色大毛片视频在线观看 | 国产精品国产三级国产av玫瑰| 国产伦精品一区二区三区四那| 91在线精品国自产拍蜜月| 永久网站在线| 亚洲精品中文字幕在线视频 | 直男gayav资源| 精品久久久久久久久av| 大香蕉久久网| 国产精品一区www在线观看| 日本午夜av视频| 久久久久久久亚洲中文字幕| 少妇裸体淫交视频免费看高清| 欧美日韩综合久久久久久| 在线免费观看不下载黄p国产| 成人国产麻豆网| av专区在线播放| 69av精品久久久久久| 精品一区二区三区人妻视频| 99久国产av精品| 亚洲18禁久久av| 欧美+日韩+精品| 99久久精品国产国产毛片| 街头女战士在线观看网站| 老女人水多毛片| 亚洲最大成人av| 亚洲成人av在线免费| 少妇熟女欧美另类| 高清日韩中文字幕在线| 夜夜看夜夜爽夜夜摸| 麻豆av噜噜一区二区三区| 国产成人免费观看mmmm| av在线亚洲专区| 亚洲欧美清纯卡通| 91aial.com中文字幕在线观看| 观看免费一级毛片| 免费看美女性在线毛片视频| 成年女人看的毛片在线观看| 一级av片app| 亚洲熟妇中文字幕五十中出| 久久精品国产亚洲av涩爱| 久久久色成人| 日本午夜av视频| 久久精品熟女亚洲av麻豆精品 | 91狼人影院| 高清欧美精品videossex| 国产黄频视频在线观看| 国产成人精品一,二区| 偷拍熟女少妇极品色| av天堂中文字幕网| 三级国产精品片| 国产精品国产三级国产av玫瑰| 国产伦在线观看视频一区| 99久国产av精品国产电影| av国产免费在线观看| 国产欧美日韩精品一区二区| 18禁在线无遮挡免费观看视频| 神马国产精品三级电影在线观看| 老师上课跳d突然被开到最大视频| 国产精品女同一区二区软件| 欧美三级亚洲精品| 91精品国产九色| 国产成人福利小说| 国产精品1区2区在线观看.| 亚洲精品成人av观看孕妇| eeuss影院久久| 成人午夜精彩视频在线观看| 成人无遮挡网站| 日本与韩国留学比较| 免费少妇av软件| 日韩av在线大香蕉| 91午夜精品亚洲一区二区三区| 亚洲精品一区蜜桃| 最近视频中文字幕2019在线8| 欧美三级亚洲精品| 亚洲国产成人一精品久久久| 国产一区二区亚洲精品在线观看| 成人亚洲欧美一区二区av| 亚洲一区高清亚洲精品| 少妇的逼好多水| 国产精品一区二区在线观看99 | 中文字幕av成人在线电影| 精品久久久久久久久久久久久| 激情 狠狠 欧美| 欧美xxⅹ黑人| 99热这里只有精品一区| 我要看日韩黄色一级片| 国产一级毛片七仙女欲春2| 欧美极品一区二区三区四区| 综合色av麻豆| 亚洲欧美精品自产自拍| 99热这里只有精品一区| 欧美性猛交╳xxx乱大交人| 国产免费福利视频在线观看| 女的被弄到高潮叫床怎么办| 淫秽高清视频在线观看| 国产午夜福利久久久久久| 精品久久久久久成人av| 国产精品爽爽va在线观看网站| 免费观看在线日韩| 亚洲久久久久久中文字幕| 欧美性猛交╳xxx乱大交人| 麻豆久久精品国产亚洲av| 99久久人妻综合| 欧美日韩视频高清一区二区三区二| 久久久久久九九精品二区国产| 天堂影院成人在线观看| 久久精品久久精品一区二区三区| 亚洲激情五月婷婷啪啪| 亚洲人成网站高清观看| 国产精品久久久久久av不卡| 黄片无遮挡物在线观看| 国产激情偷乱视频一区二区| 日本一二三区视频观看| 日韩不卡一区二区三区视频在线|