999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于KZ 濾波法的江淮地區PM2.5 濃度變化影響分析

2023-12-26 10:06:48吳明胤張付海
中國環境監測 2023年5期
關鍵詞:影響

董 昊,王 歡,吳明胤,張付海

安徽省生態環境監測中心,安徽 合肥 230071

20 世紀初,隨著我國經濟的高速發展,環境空氣污染形勢較為嚴峻,以霧霾天為代表的污染天在京津冀、長三角以及汾渭平原等區域時有發生。 細顆粒物PM2.5作為霧霾天的主要影響因子一方面會降低能見度,影響交通安全,另一方面則會通過呼吸道進入人體肺部,危害人體健康[1]。《中共中央關于制定國民經濟和社會發展第十四個五年規劃和二○三五年遠景目標的建議》要求基本消除重污染天氣,因而為進一步治理PM2.5,對已采取的管控措施效果進行評估顯得尤為重要。

國內外學者研究表明,污染源排放與PM2.5濃度的變化息息相關,是造 成PM2.5污染的根本原因,但氣象條件對大氣環境中PM2.5物理傳輸和化學轉化也有著重要影響,其中風速、風向、溫度、相對濕度以及混合層厚度等氣象因素對PM2.5濃度有著顯著影響[2-3]。 因而要量化管控措施效果,對氣象因素的扣除尤為重要。 目前,對環境空氣短期治理中管控效果評估主要采用數值模型進行定量分析,許艷玲等[4]采用WRF-CMAQ模型系統定量分析了2016 年氣象和排放因素變化導致全國PM2.5年均濃度下降幅度分別為4%和3%,王文丁等[5]使用基于嵌套網格空氣質量預報模式(NAQPMS)分析區域聯防聯控應急減排對成渝地區PM2.5濃度降低率為5%~11%。 但受限于空氣質量數值模式的性能、排放清單的不準確、化學動力學機理的復雜以及其他的不穩定因素,不同模式間的評估結果存在一定的差異性,對污染物濃度長時間序列變化趨勢分析存在局限性[6]。 因此,對長期污染物濃度變化趨勢分析采用Kolmogorov-Zurbenko(KZ)濾波法被更廣泛使用,其原理是利用污染物濃度變化的頻域特性對其時間序列進行不同時間尺度分離,并通過逐步回歸的方法得到氣象條件對污染物濃度變化的影響。

安徽省位于長江三角洲地區,有八百里的沿江城市群和皖江經濟帶,尤其是江淮之間及沿江區域近年來經濟高速發展,城市化、工業化程度進展較快。 本文以沿江地區的合肥、蕪湖和馬鞍山3 個城市為研究對象,分析2018—2020 年各城市PM2.5和氣象參數日數據,從而量化“十三五”期間氣象條件對沿江地區的PM2.5濃度長期變化的影響,以期有助于探討經濟高速發展下城市減排措施的實施效果。

1 研究方法

1.1 觀測數據

本研究選取2018 年1 月1 日—2020 年12 月31 日合肥、蕪湖和馬鞍山3 個城市共計10 個國控站點的監測數據(具體監測站點信息見表1),根據《環境空氣質量評價技術規范(試行)》(HJ 663—2013)[7]和《環境空氣質量標準》 (GB 3095—2012)[8]評價各城市PM2.5日均濃度,進行KZ 濾波統計分析。 同時,為保證PM2.5監測數據與氣象數據相匹配,選取同站點監測的氣象參數,包括日平均氣壓(P)、日平均氣溫(T)、日平均相對濕度(RH) 和日平均風速(v),并進行數據整理。

表1 空氣質量監測站點信息及分布Table 1 Information and distribution of air quality monitoring stations

1.2 研究方法

濾波統計是一種時間序列分析方法,其原理是利用原始參數變化的頻域特性對其時間序列進行不同時間尺度的分離。 對某一原始參數的時間序列A(t),其模型可表示為

式中:W(t)是參數短期分量;S(t)是參數季節分量;e(t)是參數長期分量。 其中,針對大氣環境中PM2.5濃度的時間序列變化,其短期分量主要受短時天氣系統影響,而排放源除了在特定節假日或其他特殊事件會顯著影響濃度變化外,其余時段變化總體較小。 季節分量則是地球公轉導致太陽角度的變化,使得氣象條件的影響出現季節性變化,例如降水、季風以及氣溫的變化,同時因季節變化導致的污染排放也是影響之一。 長期分量變化的影響因素較多,主要包括區域內氣候變化趨勢、污染物的排放總量、經濟發展狀況以及當地相關政策[9]。

KZ 濾波是一種低通濾波,是m 個點的滑動平均濾波器的p 次迭代。 其計算公式:

式中:Yi為經KZ 濾波后的時間序列;m 為滑動窗口長度(m=2k+1);i 為序列的時間間隔;j 為滑動窗口變量;k 為Ai在濾波時其兩端的滑動窗口長度。 將濾波后的時間序列作為下一次濾波的原始時間序列進行計算,根據計算需求,重復執行p次,得到最終的濾波結果用KZ(m,p)表示。 不同尺度的濾波是通過調整參數m 和p 來控制的,其中參數m 和p 與有效濾波寬度N 滿足如下公式[10]:

根據RAO 等[11]和吳宜航等[12]的研究結果,選取KZ(365,3)和KZ(15,5)進行濾波處理,有效濾波分別為33 d 和632 d。 其中,通過KZ(365,3)濾波可以提取時間序列的長期分量,而經過KZ(15,5)濾波后的時間序列可定義為基線分量XB(t),代表季節分量和長期分量的總和,具體公式:

根據式(1)、式(4)和式(5)可以得到季節分量和短期分量的計算公式:

為進一步消除氣象因素對PM2.5濃度的影響,本研究通過以氣象參數為自變量的多元線性回歸,分別就短期分量和基線分量構建氣象要素與PM2.5的關系模型,并計算各分量回歸模型的解釋方差:

式中:VE 代表解釋方差;varA(t)表示污染物原始序列方差;varε(t)表示殘差序列方差。 解釋方差越大,氣象要素對原始序列的影響越大,解釋能力越強。

各分量與回歸模型預測值的殘差序列即為剔除了氣象參數影響的PM2.5濃度序列,即由污染排放導致的PM2.5濃度變化。 其中,原始序列的總殘差序列ε(t)是短期分量回歸殘差εM(t) 與基線分量回歸殘差εB(t) 之和,計算公式:

繼續對總殘差序列做KZ(365,3)濾波處理,消除其他未參與擬合的氣象因子的影響,獲得受排放變化影響的PM2.5濃度長期變化分量εL(t) 。將變化分量與原始序列的長期分量均值相疊加,得到僅由污染源排放影響的PM2.5濃度長期分量時間序列eadj(t):

基于排放源和氣象條件對PM2.5濃度影響的區分,計算2018—2020 年合肥、蕪湖和馬鞍山扣除氣象條件后PM2.5濃度,對比分析濃度變化,從而評估分析時段氣象條件和污染源減排措施對PM2.5濃度趨勢變化的貢獻。

2 結果與討論

2.1 KZ 濾波后PM 2.5 的時間序列特征

根據《環境空氣質量標準》(GB 3095—2012)和《環境空氣質量指數(AQI)技術規定(試行)》(HJ 633—2012)中的規定,PM2.5的日均二級質量濃度限值為 75 μg/m3, 日均濃度超過150 μg/m3即為重污染。 圖1 為PM2.5日濃度原始序列,由圖1 可知:合肥、蕪湖和馬鞍山PM2.5的日變化趨勢較為一致,日超標天主要集中在冬季,其中PM2.5日均濃度超過150 μg/m3的重污染天占總天數的比例分別為0.6%、1.1% 和0.7%。 統計分析2018—2020 年3 個城市的重污染天數變化情況,發現3 個城市重污染天數均呈現逐年減少趨勢;對比城市間日均變化,合肥和蕪湖PM2.5的日超標天數相對較多,超標率分別為10.3%和10.6%,馬鞍山PM2.5污染程度相對較輕,日超標率為8.6%。

圖1 合肥、蕪湖和馬鞍山PM 2.5 的原始時間序列A(t)Fig.1 The original time series A(t) of PM 2.5 in Hefei,Wuhu and Ma'anshan City

圖2 ~圖4 為合肥、蕪湖和馬鞍山3 個城市代表站點的PM2.5日濃度經KZ 濾波法處理后得到的時間序列。 短期分量變化如圖2 所示,合肥、蕪湖和馬鞍山的高頻變化較為一致,振幅范圍分別為-65 ~+101、-76 ~+135 和-76 ~+122,反映了周期小于33 d 的PM2.5濃度波動情況,考慮污染排放量在短期內總體較為穩定,因而短期分量的波動主要是受本地天氣變化影響。 同時,3 個城市在冬春季(11 月至次年3 月)波動明顯強于其余時段,說明不同季節的天氣影響存在差異,冬春季受本地天氣變化影響更為顯著。 因而在冬春季進行人影等氣象作業,能更易達到本地PM2.5濃度削減的目的。

圖2 合肥、蕪湖和馬鞍山PM 2.5 經KZ 濾波分解后的短期分量序列W(t)Fig.2 The short-term component W(t) of PM 2.5 by the KZ filter in Hefei,Wuhu and Ma'anshan City

去除短期分量的影響后,季節分量(圖3)反映的是周期為33—632 d 的PM2.5濃度波動,主要是受本地排放源生產周期和氣象季節條件變化影響,合肥、蕪湖和馬鞍山的季節分量變化波動基本一致,呈現明顯周期季節特征。 其中,各城市濃度峰值主要集中在12 月至次年1 月,主要原因一方面是氣溫偏低,大氣混合層較低,相對濕度較高,天氣系統以靜穩為主,大氣水平和垂直擴散能力差,有利于本地污染物的累積[13],另一方面江淮之間的冬季主導風向為偏北風,京津冀及周邊地區冬季采暖導致的區域性污染輸送影響較為顯著[14]。 但在春季PM2.5濃度也會出現短時的上升,如2018 年5 月下旬,這是由于上游沙源地起沙后對江淮之間區域的輸送影響[15],其中受地理位置的影響,合肥受沙塵影響更為明顯,PM2.5濃度的波動較蕪湖和馬鞍山更為顯著。 而在7—8月,受梅雨季影響疊加夏季有利的氣象條件,3 市降水多、溫度較高,大氣垂直對流較為活躍,有利于污染物的擴散和沉降,故均在此時段出現PM2.5濃度波谷。 此外,不同城市間的波動周期也存在一定的差異性,如合肥和蕪湖在2018 年11 月至2019 年1 月的峰值出現分叉,而馬鞍山未出現;波谷出現先后順序為馬鞍山>蕪湖>合肥。

圖3 合肥、蕪湖和馬鞍山PM 2.5 經KZ濾波分解后的季節分量序列S(t)Fig.3 The seasonal component S(t) of PM 2.5 by the KZ filter in Hefei,Wuhu and Ma'anshan City

圖4 為濾波后的長期分量序列e(t),其濾波周期是632 d,其濃度波動反映的是區域污染物排放總量和氣候變化等因素影響結果。 各城市的變化趨勢總體較為一致,呈現下降趨勢,其中合肥和蕪湖在2018 年11 月至2019 年3 月出現輕微波動,隨后下降趨勢趨于明顯,而這種變化趨勢是無法直接從原始序列中解讀的,進一步說明采用時間序列模型進行分析的必要性[16]。 合肥、蕪湖和馬鞍山長期分量下降速率分別為3.6、3.8 和2.6 μg/(m3·a),PM2.5年均值下降速率為3.3、4.3 和2.3 μg/(m3·a),對比二者發現,氣象因素對合肥和蕪湖PM2.5濃度下降呈現抑制作用,而對馬鞍山則是促進作用。 此外,2018—2019 年合肥和蕪湖PM2.5的長期分量明顯高于馬鞍山,這一方面是合肥和蕪湖的本地污染源排放更為嚴重,包括工業企業廢氣、機動車尾氣以及建筑工地揚塵等[17];另一方面江淮之間區域在氣候上總體會受到北方PM2.5污染傳輸影響,相對于馬鞍山,合肥和蕪湖的地理位置更為偏北,受PM2.5污染傳輸影響更為顯著。 但隨著合肥和蕪湖的工業結構的轉型,電廠等燃煤企業的超低排放,各市間的污染排放量逐漸接近,PM2.5的濃度水平在2020 年已經與馬鞍山基本趨于一致。

圖4 合肥、蕪湖和馬鞍山PM 2.5 經KZ濾波分解后的長期分量序列e(t)Fig.4 The long trend com ponent e(t) of PM 2.5 by the KZ filter in Hefei,Wuhu and Ma'anshan City

為進一步驗證KZ 濾波的準確性,選取各分量的協方差之對原始序列方差的貢獻來檢驗時間序列分解效果,貢獻越小,說明KZ 濾波分解效果越好,各分量越獨立。 由表2 可知,3 個城市的原始序列經KZ 濾波后,協方差對總方差的貢獻均低于10%,說明各城市分解后的分量之間相對較為獨立,相互間的干擾較小。 而分析各分量方差對總方差的貢獻情況,可以發現各市中短期分量的貢獻最大,均超過50%,其次是季節分量,貢獻率范圍為40.4%~41.6%,而長期分量貢獻占比最少,其中,馬鞍山的貢獻率僅為0.9%,這說明PM2.5濃度的變化主要受污染源排放和氣象條件的短期和季節性變化影響,這與秦人潔等對石家莊、保定和張家口PM2.5的分析結果較為一致[18]。

表2 合肥、蕪湖和馬鞍山PM 2.5 各分量方差對總方差的貢獻率Table 2 Contribution of each component to the total variance of PM 2.5 in Hefei,Wuhu and Ma'anshan City

2.2 氣象因素對PM 2.5 不同分量的影響

多項研究表明,PM2.5的濃度變化受氣象因素影響較為顯著,為進一步探究經KZ 濾波后不同分量下PM2.5與各氣象因子的相關性,本研究使用Spearman 相關系數法進行分析,結果見表3。相對于長期分量中,PM2.5與氣壓呈現弱負相關,在季節分量中則為顯著正相關,這與王嫣然等在北京地區對不同季節PM2.5與氣象因素的分析結果一致[19];與氣溫的相關系數在季節和長期分量中呈現負值,而在短期分量中為正值,這主要是氣溫上升會促進大氣垂直運動,有利于PM2.5的垂直擴散和稀釋,但若在短期內有逆溫層存在,隨著氣溫的升高,有利于大氣穩定,大氣混合層高度偏低,從而使得近地面污染物聚集[20];相對濕度和風速與PM2.5有著顯著相關性[21],其中相對濕度在長期分量中與PM2.5呈現負相關,但在合肥的短期和季節分量出現正相關,這有可能是相對濕度較高導致PM2.5吸濕增長;風速對PM2.5的影響主要是通過水平擴散,隨著風速增大,PM2.5的傳輸速率增加,污染物濃度下降,因而風速在短期和季節分量上與PM2.5呈現負相關,但城市的基礎設施建設作為顆粒物的重要排放源,在風速增大的時候,會造成建筑工地揚塵擴散,使得空氣中顆粒物濃度上升,這可能是風速與PM2.5長期分量呈現正相關的原因[20]。 綜上,各氣象因子在不同時間尺度下對PM2.5的影響作用偏差較大,應綜合考慮不同時間尺度下氣象條件對污染物的影響。

表3 合肥、蕪湖和馬鞍山PM 2.5 與氣象因子的Spearm an 相關系數Table 3 Spearm an's correlation coefficients between PM 2.5 and meteorological factors in Hefei,Wuhu and Ma'anshan City

2.3 基于氣象因素調整PM 2.5 的濃度變化

為進一步消除氣象因素對PM2.5濃度變化的影響,根據上述氣象因子對經KZ 濾波后PM2.5時間序列進行多元線性回歸分析,并構建統計模型。由表4 可知,各城市均呈現出長期分量的解釋方差最高、季節分量次之、短期分量最低的變化趨勢。其中,各城市的長期分量和季節分量解釋方差為87.9%~99.2%,可以認為長期分量和季節分量基本能通過選取的4 個氣象因子變量來解釋[20]。

表4 合肥、蕪湖和馬鞍山PM 2.5 短期分量、季節分量和長期分量的回歸模型分析Table 4 Regression model of PM 2.5 concerning the short-term component,seasonal component,and long-trend component by the stepw ise regression in Hefei,Wuhu and Ma'anshan City

圖5 是依據回歸模型和殘差分析后繪制的PM2.5濃度長期分量變化趨勢圖。 由圖5 可見:各城市消除氣象因素影響后長期分量總體有所下降,尤其是2018—2019 年,說明氣象條件整體不利于PM2.5的改善,這與梅梅等在京津冀采用城市大氣污染數值預報系統進行分析得到2019 年氣象條件對PM2.5濃度上升貢獻9.96%的結果較為一致[22];而在2020 年下半年,氣象條件則對PM2.5的影響呈現改善作用。

圖5 合肥、蕪湖和馬鞍山經氣象修正前后PM 2.5 濃度長期分量變化趨勢Fig.5 Variations of long-trend component of PM 2.5 concentration before and after meteorological adjustm ent in Hefei,Wuhu and Ma'anshan City

對比各城市2018 年和2020 年長期分量修正前后的PM2.5濃度均值,其中2018 年與2020 年長期分量修正前PM2.5濃度均值的差值即為污染源與氣象條件的總貢獻,而長期分量修正后的差值為污染源的總貢獻,具體結果見表5。 同2018年相比,2020 年合肥、蕪湖和馬鞍山PM2.5污染均顯著改善,污染源減排措施仍是PM2.5污染改善的主要因素,其中氣象條件對合肥、蕪湖和馬鞍山PM2.5改善的貢獻率分別為 1.0%、 7.2% 和17.3%。

表5 合肥、蕪湖和馬鞍山減排措施和氣象條件對PM 2.5 濃度改善的貢獻情況Table 5 Contribution of emission reduction measures and meteorological conditions to the improvement of PM 2.5 in Hefei,Wuhu and Ma'anshan City

本研究基于KZ 濾波法評估了氣象和污染源對3 個城市PM2.5的影響情況,但仍有其不確定性。 主要包括:①監測站點對城市的代表性有限,不同區域的PM2.5監測數據和氣象條件都存在差異,尤其是近地面氣象參數,受各種條件的制約,本研究僅選取10 個站點數據進行了均值處理;②基于多元線性模型進行擬合時,采用的擬合因子不夠全面,擬合方程僅考慮氣象因子與PM2.5濃度之間的關系,未考慮大氣環境中二次生成對PM2.5的貢獻。 ③貢獻分析采用的是長期分量的計算結果,對于短期和季節分量中氣象條件影響未進行剝離。

2.4 基于數值模式的情景分析

嵌套網格空氣質量模式系統(NAQPMS)模式采用中尺度天氣預報模式(WRF)輸出的氣象要素場作為的動力驅動,美國國家環境預報中心的NCEP 再分析數據為初始條件和邊界條件,基于清華大學MEIC 清單為污染排放源清單,設置模擬區域為三重區域嵌套,網格分辨率分別為27、9、3 km,第三層網格覆蓋安徽省全境。 本研究采用NAQPMS 模型進行3 個城市的情景模擬,保持排放源清單不變,第1 組模擬方案采用2018 年的氣象場,第2 組模擬方案采用2020 年的氣象場,兩者模擬PM2.5濃度值為氣象條件帶來的濃度變化。

采用標準化分數偏差(MFB)對NAQPMS 模型模擬結果進行評估[式(11)]。

式中:Cm為模型模擬值;Co為實際監測值。

數字模型模擬結果顯示,2018 年合肥、蕪湖和馬鞍山的PM2.5模擬值與實測值的平均偏差分別為11.9,13.6 和7.8,MFB 分別為23.1%,14.9%和28.2%,模擬結果在理想范圍內[23],而造成模擬值偏高的原因可能是源清單的不確定性和模型本身化學機理的不完善。

為縮小模式模擬的差異,根據2018 年數據對3 個城市采用最小二乘法,構建PM2.5濃度模擬值與實測值的線性擬合關系(圖6)。 基于擬合方程對2020 年3 個城市模擬值進行調整,調整后的PM2.5濃度為排放強度與2018 年一致污染物能達到的濃度水平。

圖6 2018 年合肥、蕪湖和馬鞍山PM 2.5 濃度模擬值和實測值線性擬合Fig.6 Linear regression of simulated and observed PM 2.5 concentrations in Hefei,Wuhu and Ma'anshan City in 2018

圖7 為NAQPMS 模式模擬的氣象和污染源變化對2020 年3 個城市PM2.5濃度的影響。 結果顯示,相比于2018 年,2020 年合肥和馬鞍山PM2.5濃度下降的主要貢獻來自于減排措施,與KZ 濾波的結論較為一致;蕪湖PM2.5濃度下降貢獻中,減排措施和氣象條件分別占比42.5%和57.5%,較KZ 濾波統計算法結果有所偏差,其原因可能是實際減排控制措施的執行效果好于減排預期效果,致使數值模型對減排措施存在低估。

圖7 基于NAQPMS 模式合肥、蕪湖和馬鞍山減排措施和氣象條件對PM 2.5 濃度改善的貢獻情況Fig.7 Contribution of emission reduction measures and meteorological conditions to the improvement of PM 2.5 in Hefei,W uhu and M a'anshan City based on NAQPM S

3 結論

1)江淮之間3 個城市PM2.5濃度各分量對原始序列總方差的貢獻和均超過90%,KZ 濾波的分離尺度較為合理,能較為科學構建區域內PM2.5濃度時間序列模型,其中各分量中,短期分量對總方差的貢獻最大,季節分量次之,長期分量貢獻最少,說明PM2.5濃度時間序列的波動主要是受污染源和氣象條件的短期和季節影響。

2)PM2.5濃度與氣壓、溫度、相對濕度和風速均呈現一定相關性。 其中,在不同時間尺度下,氣象因子對PM2.5的影響存在顯著差異。

3)依據回歸模型和殘差分析后,修正氣象條件對長期分量的影響,KZ 濾波統計結果顯示2018—2020 年氣象條件對江淮之間區域PM2.5污染改善影響存在波動,同2018 年相比,2020 年江淮之間PM2.5污染均顯著改善,減排措施的貢獻率為82.7%~99.0%;數值模式情景分析結果顯示減排措施的貢獻率為42.5% ~104.8%;因而,污染源減排措施是江淮之間PM2.5污染改善的主要因素。

猜你喜歡
影響
是什么影響了滑動摩擦力的大小
哪些顧慮影響擔當?
當代陜西(2021年2期)2021-03-29 07:41:24
影響大師
沒錯,痛經有時也會影響懷孕
媽媽寶寶(2017年3期)2017-02-21 01:22:28
擴鏈劑聯用對PETG擴鏈反應與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
基于Simulink的跟蹤干擾對跳頻通信的影響
如何影響他人
APRIL siRNA對SW480裸鼠移植瘤的影響
對你有重要影響的人
主站蜘蛛池模板: 99国产精品一区二区| 久久久久国产一级毛片高清板| 在线免费a视频| 在线看国产精品| 91蜜芽尤物福利在线观看| 蜜桃视频一区| 香蕉99国内自产自拍视频| 亚洲成年人片| 亚洲精品无码在线播放网站| 欧美日韩亚洲国产| 亚洲国产理论片在线播放| 亚洲av无码久久无遮挡| 日本尹人综合香蕉在线观看| 国产午夜人做人免费视频| 2020久久国产综合精品swag| 欧美成人免费午夜全| 欧美国产日产一区二区| 国产杨幂丝袜av在线播放| 日本午夜精品一本在线观看| 亚洲制服丝袜第一页| 亚洲日产2021三区在线| 国产成人1024精品下载| AV熟女乱| 国产成人高清精品免费软件 | 日本在线视频免费| 亚洲一区第一页| 国产AV无码专区亚洲精品网站| 女人18毛片久久| 国产精品一区二区不卡的视频| 国产精品私拍在线爆乳| 蜜臀AVWWW国产天堂| 亚洲三级片在线看| 在线观看亚洲人成网站| 国产欧美一区二区三区视频在线观看| 男人的天堂久久精品激情| 色综合日本| 98超碰在线观看| 成人久久18免费网站| 国产乱子伦手机在线| 热久久这里是精品6免费观看| 最新国产成人剧情在线播放 | 国产永久在线视频| 99九九成人免费视频精品| 久久无码av一区二区三区| 少妇露出福利视频| 丁香五月激情图片| 国产精品三区四区| 精品国产欧美精品v| 国产精品999在线| 国产亚洲视频在线观看| 国产91全国探花系列在线播放| 国产h视频免费观看| 毛片在线播放网址| 凹凸国产熟女精品视频| 被公侵犯人妻少妇一区二区三区| 五月天福利视频 | 国产黄网站在线观看| 日本道中文字幕久久一区| 人妻精品久久无码区| 国产系列在线| 免费国产好深啊好涨好硬视频| 四虎永久免费地址在线网站| 色吊丝av中文字幕| 午夜丁香婷婷| 欧美在线三级| 4虎影视国产在线观看精品| 国语少妇高潮| 欧美精品一二三区| 欧美日韩综合网| 国产第一页第二页| 蜜臀AVWWW国产天堂| 久久精品亚洲专区| 国产精品免费p区| 亚洲v日韩v欧美在线观看| 亚洲啪啪网| 欧美翘臀一区二区三区| 在线欧美国产| 欧美黑人欧美精品刺激| 国产手机在线观看| 国产成人亚洲综合a∨婷婷| 国产精品久久久久婷婷五月| 91综合色区亚洲熟妇p|