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

中國霧霾污染的空間計量分析

2018-09-21 05:42:58林楚海
統(tǒng)計與決策 2018年16期
關(guān)鍵詞:分析模型

林楚海

(暨南大學(xué) 經(jīng)濟(jì)學(xué)院,廣州 510000)

0 引言

霧霾,即霧和霾的組合詞。霧是由大量懸浮在近地面空氣中的微小水滴或冰晶組成的氣溶膠系統(tǒng),而霾是指由于灰塵、硫酸、硝酸、有機碳?xì)浠衔锏攘W討腋≡诳諝庵校顾侥芤姸刃∮?000米一種現(xiàn)象。霧霾天氣主要由空氣中的二氧化硫、氮氧化物和可吸入顆粒物這三項引起,它們與霧氣結(jié)合在一起,會讓天空瞬間變得陰沉灰暗,也對人體和社會帶來各方面的影響。2014年,中國31個省份中,PM10年平均濃度達(dá)到國家一級標(biāo)準(zhǔn)的只有海南省,而SO2達(dá)到國家一級標(biāo)準(zhǔn)的也只有7個省份,從整體上看,霧霾污染依然不容樂觀。再加上霧霾的源頭多種多樣,比如汽車尾氣、工業(yè)排放、建筑揚塵、垃圾焚燒等,因此,如何在眾多因素中快速找到造成霧霾的主導(dǎo)因素并采取有效的治理措施,成為中國政府和居民急需解決的難題。基于此,本文將選取中國31個省份2014年三種霧霾主要組分——PM10、SO2、NO2的年平均濃度作為目標(biāo)變量,并利用LASSO高維變量選擇方法從22個社會影響因子中篩選出非零變量進(jìn)行空間回歸分析,從而揭示中國內(nèi)地霧霾污染的真正來源和空間影響因素。

1 研究方法

1.1 空間計量分析

1.1.1 空間權(quán)重矩陣

空間權(quán)重矩陣即地理單元間的空間結(jié)構(gòu)關(guān)系的一種表達(dá)形式,描述了各單元間的鄰近關(guān)系。空間權(quán)重矩陣W一般表示為一個n維矩陣,其基本形式可表示如下:

其中,wij表示區(qū)域i和j的鄰近關(guān)系,一般分為基于鄰近概念的空間權(quán)重矩陣和基于距離的空間權(quán)重矩陣。其定義如下:

1.1.2 全局空間自相關(guān)

全局空間自相關(guān)描述的是研究范圍內(nèi)所有單元之間的整體空間關(guān)系。一般采用Moran’s I統(tǒng)計量來描述整個研究區(qū)域上所有空間單元之間的平均關(guān)聯(lián)程度及其顯著性。其計算公式如下:

其中,xi表示地區(qū)i的觀測值,,n為地區(qū)數(shù),wij為空間權(quán)重矩陣。Moran’s I統(tǒng)計量的取值范圍為[-1,1]。I越接近1,表示區(qū)域間空間正相關(guān)的程度越強;當(dāng)I越接近-1時,表示區(qū)域間空間負(fù)相關(guān)的程度越強;當(dāng)I接近0時,則表示區(qū)域間不存在空間自相關(guān)性。

1.1.3 局域空間自相關(guān)

局域空間自相關(guān)統(tǒng)計量就是用來識別隨空間位置不同而存在的不同空間關(guān)聯(lián)模式,從而觀察空間局部不平穩(wěn)性。本文的局域空間自相關(guān)分析將采用局部Moran’s I指數(shù)繪制的集聚地圖和顯著性地圖。

局部Moran’s I指數(shù)的計算公式如下:

其中,Ii表示i地區(qū)與其相鄰區(qū)域之間的空間相關(guān)程度,其他符號的含義同式(4),這里不再贅述。當(dāng)Ii>0時,表示第i個地區(qū)某種屬性值與其周圍地區(qū)的屬性值呈空間正相關(guān),表現(xiàn)為高高集聚或低低集聚;當(dāng)Ii<0時表示負(fù)相關(guān),表現(xiàn)為高低集聚或低高集聚。

1.1.4 空間回歸分析

(1)空間滯后模型:假定空間變量的空間依賴性僅通過因變量產(chǎn)生,即某一空間對象上的因變量不僅與同一對象上的自變量有關(guān),還與相鄰對象的因變量有關(guān),其形式為:

其中,y為n×1階截面數(shù)據(jù)被解釋變量,W為空間權(quán)重矩陣,Wy即為被解釋變量的空間滯后向量,其系數(shù)ρ為空間自回歸系數(shù),反映了相鄰區(qū)域?qū)Ρ緟^(qū)域的作用,當(dāng)ρ顯著時,說明因變量間存在明顯的空間依賴性,當(dāng)ρ>0時,表明相鄰空間區(qū)域的空間相互作用表現(xiàn)為“溢出效應(yīng)”;反之,則表現(xiàn)為“負(fù)效應(yīng)”或“擴(kuò)散效應(yīng)”。X即為n×k外生解釋變量矩陣,其中k為解釋變量的個數(shù)。ε為n×1階隨機誤差向量。

(2)空間誤差模型:在某些情況下,空間依賴可能是某些對局域內(nèi)多個單元均有影響的遺漏變量造成的,從而造成區(qū)域“一損俱損,一榮俱榮”的效果,這時候需要考慮誤差項的空間自回歸過程。其形式為:

其中,λ是誤差項的空間自回歸系數(shù),Wε是誤差項的空間滯后向量,μ為不相關(guān)的、均值為0且同方差的誤差項。若λ顯著,則說明了鄰近區(qū)域的關(guān)于因變量的誤差沖擊對本地區(qū)的觀測值具有顯著性影響。

(3)空間通用模型:空間通用模型綜合了空間滯后模型和空間誤差模型,其形式為:

其中,空間權(quán)重矩陣W1與W2可以相同,也可以不同,其他符號含義與上述兩個模型相同,這里不再贅述。

(4)空間相關(guān)性檢驗與模型選擇:一般可通過兩個拉格朗日乘數(shù)LMERR、LMLAG及其穩(wěn)健的R-LMERR、R-LMLAG等統(tǒng)計量來實現(xiàn)。當(dāng)LMLAG比LMERR顯著,且R-LMLAG顯著而R-LMERR不顯著時,則選擇空間滯后模型更為合適;相反地,如果LMERR比LMLAG顯著,且R-LMERR顯著而R-LMLAG不顯著時,則選擇空間誤差模型更為合適。其他可供參考的模型比較準(zhǔn)則還有對數(shù)似然值、赤池信息準(zhǔn)則(AIC)和施瓦茨準(zhǔn)則(SC)等。當(dāng)對數(shù)似然值越大,AIC和SC值越小,則模型擬合效果越好。

1.2 LASSO高維變量選擇

設(shè)有p個自變量 x1,x2,…,xp和因變量y,并建立如下的線性回歸模型:

其中,α為常數(shù)項,β1,β2,…,βp為自變量的回歸系數(shù),ε為隨機擾動項。設(shè)(xi1,xi2,…,xip;yi),i=1,2,…,n 為n組觀測值,假定數(shù)據(jù)已經(jīng)經(jīng)過中心化和標(biāo)準(zhǔn)化,即,…,βp)T,則 β 的最小二乘估計為:

1996年,Tibshirani提出了LASSO變量選擇方法,其基本思想是在最小二乘的基礎(chǔ)上施加l1懲罰項,將一些變量的系數(shù)恰好壓縮為零,從而實現(xiàn)變量選擇。其表達(dá)式如下:

其中,s≥0是一個懲罰參數(shù),控制著壓縮的程度,s越小,壓縮程度越強,會有更多的系數(shù)被壓縮至零。記最小二乘估計得到的參數(shù)估計為實現(xiàn)壓縮,即有些系數(shù)會變?yōu)?,這些等于0的變量就會被刪除,從而達(dá)到變量選擇的目的。其等價形式如下:

若設(shè)計陣正交,即XTX=I,LASSO的參數(shù)估計為:

2 實證分析

2.1 數(shù)據(jù)說明

本文選取中國31個省份2014年的PM10、SO2、NO2的年平均濃度(微克/立方米)作為被解釋變量,同時從產(chǎn)業(yè)結(jié)構(gòu)、人口因素、固定資產(chǎn)投資和房地產(chǎn)、能源因素、城市建設(shè)五個大方面選擇22個指標(biāo)作為自變量,數(shù)據(jù)均來源于各省2015年的《統(tǒng)計年鑒》和《環(huán)境統(tǒng)計公報》,如下頁表1所示。

采用ArcGis 10.0軟件計算空間權(quán)重矩陣,相鄰規(guī)則采用Rook相鄰。另外,由于海南島遠(yuǎn)離陸地,與其他省份沒有邊界相鄰,為了避免出現(xiàn)“孤島”現(xiàn)象,將鄰居個數(shù)最低值設(shè)為1,從而使得每個省份至少都有一個鄰居相鄰。

表1 指標(biāo)含義與說明

2.2 全局空間自相關(guān)分析

通過Geoda軟件進(jìn)行全局空間自相關(guān)分析,得到31個省份2014年的PM10、SO2、NO2年平均濃度對應(yīng)的全局Moran’s I指數(shù),如表2所示。結(jié)果表明,PM10、SO2、NO2年平均濃度對應(yīng)的全局Moran’s I指數(shù)分別為0.5377、0.3501、0.5044,并且對應(yīng)的P值均小于0.05,說明31個省份之間的PM10、SO2、NO2的空間自相關(guān)性非常顯著,且呈正的空間自相關(guān),而非完全隨機分布。

表2 全局Moran’s I指數(shù)及相關(guān)統(tǒng)計指標(biāo)

2.3 局域空間自相關(guān)分析

在現(xiàn)實中,各地區(qū)霧霾組分受到的空間影響的模式完全一致一般是很少見的,為了進(jìn)一步探究各省份與相鄰省份之間的空間影響關(guān)系,本文繪制了各省份PM10、SO2、NO2的顯著性地圖和聚類地圖(圖略)。在0.05的顯著性水平下,華北地區(qū)的北京、天津、山西、河北以及華東地區(qū)的山東、華中地區(qū)的河南這六個省份的PM10年平均濃度具有顯著性的空間正自相關(guān),呈“高高集聚”的分布形態(tài);而華南地區(qū)的廣東、廣西和西南地區(qū)的云南、貴州的PM10年平均濃度也具有顯著性的正自相關(guān),但呈“低低集聚”的分布形態(tài)。其他省份的PM10年平均濃度的空間自相關(guān)性則不顯著,呈隨機分布的形態(tài)。可以發(fā)現(xiàn),中國PM10的高濃度聚集區(qū)主要分布在華北地區(qū)以及華東、華中地區(qū)的個別省份,這幾個地區(qū)是中國人口較為密集和工業(yè)生產(chǎn)較為頻繁的地區(qū),工業(yè)氣體廢棄物的排放和眾多人口冬季的采暖燃煤可能是造成這幾個地區(qū)PM10濃度高的主要原因。

對于SO2來說,山東、河北、河南這三個省份的SO2年平均濃度具有顯著性的空間正相關(guān)關(guān)系,呈現(xiàn)“高高集聚”的空間布局。而北京和內(nèi)蒙古雖然都與高SO2的河北相鄰,但卻均呈現(xiàn)出“低高集聚”,對于北京來說,這可能與政府近幾年的污染治理息息相關(guān),SO2是與燃煤直接相關(guān)的大氣污染物。自1998年以來,在燃煤鍋爐清潔能源改造上,北京城六區(qū)累計完成約5.13萬蒸噸燃煤鍋爐改造,基本實現(xiàn)了城六區(qū)無燃煤鍋爐,因此,北京的SO2年平均濃度才能在高SO2的相鄰區(qū)域下保持較低的濃度,也從側(cè)面反映出工業(yè)污染治理對環(huán)境帶來的巨大效應(yīng)。而對于內(nèi)蒙古來說,其SO2年平均濃度之所以能保持較低的濃度,可能與其地理優(yōu)勢有關(guān),內(nèi)蒙古北有大青山、東南部有蠻漢山、西南及南部則為平原,且由東北向西南逐漸傾斜,因此,形成了阻礙SO2擴(kuò)散的低氣壓帶,對河北的SO2的空間傳輸形成了較強的阻礙作用。而廣東屬于沿海省份,海陸風(fēng)對霧霾組分的擴(kuò)散和稀釋起了很大的作用,再加上自身污染沒有北方省份嚴(yán)重,因此,其SO2年平均濃度依然呈現(xiàn)“低低集聚”,即與周圍省份相似,SO2年平均濃度都相對較低。

另外,北京、天津、河北三省的NO2年平均濃度呈現(xiàn)“高高集聚”的分布形態(tài),說明這三個省份也是NO2的高濃度地區(qū),而廣東的NO2年平均濃度則相對較低,呈“低低集聚”分布。

2.4 基于LASSO方法的變量選擇

由全局空間自相關(guān)分析和局域空間自相關(guān)分析可以發(fā)現(xiàn),中國各省份之間的霧霾組分確實存在空間效應(yīng),因此,需要對其進(jìn)行空間回歸分析,以探討哪種回歸模型更適合于三種霧霾組分。在建模之前,先使用LASSO方法對變量進(jìn)行篩選,以使得模型更加簡潔和更有解釋性。為了消除各變量量綱的影響,先對各變量取對數(shù),其中,X1(第二產(chǎn)業(yè)增加值占GDP比重)、X3(城鎮(zhèn)人口比重)、X11(建成區(qū)綠化覆蓋率)、X15(生活垃圾無害化處理率)、X16(森林覆蓋率)由于是百分比的形式,因此不進(jìn)行對數(shù)化處理。然后,對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化、中心化處理。最后通過LASSO方法得到各組分的非零變量如下頁表3所示。

2.5 空間回歸分析

接下來,利用R語言的spdep安裝包,根據(jù)LASSO篩選出的非零變量對三種霧霾組分分別進(jìn)行空間回歸分析,以探討最優(yōu)的擬合模型。

2.5.1 PM10的空間回歸分析與檢驗

建立PM10的回歸模型,并進(jìn)行診斷性檢驗,結(jié)果如表4所示。結(jié)果表明,在lnPM10的普通最小二乘估計中,X1、lnX10、X16的回歸系數(shù)對應(yīng)的P值均小于0.05,均通過顯著性檢驗。從拉格朗日乘數(shù)檢驗結(jié)果來看,LMLAG對應(yīng)的P值為0.00214,遠(yuǎn)小于LMERR對應(yīng)的P值0.12300,而且R-LMLAR的P值為0.00194<0.05,通過了顯著性檢驗,而R-LMERR的P值為0.10980,在5%的顯著性水平下不顯著。因此,應(yīng)該選擇空間滯后模型較為合適。從空間滯后模型的擬合效果來看,擬合優(yōu)度為74.669%,比最小二乘估計的擬合優(yōu)度(60.470%)高,而且其對數(shù)似然值(-23.75857)比最小二乘估計的(-29.09351)大,AIC和SC也比最小二乘估計的小,因此,采用空間滯后模型是最優(yōu)的,其表達(dá)式如下:

表3 LASSO方法篩選出的非零變量

從模型的回歸系數(shù)來看,各變量依然對PM10具有顯著性的影響,并且其空間回歸系數(shù)ρ=0.582,對應(yīng)的P值為0.00001,也具有很強的顯著性,說明中國各省份之間的PM10濃度具有很強的空間溢出效應(yīng),即當(dāng)空間滯后值每增加1%時,本地的PM10濃度將平均上升0.582%。從社會影響因子來看,當(dāng)?shù)诙a(chǎn)業(yè)占GDP比重每上升1%時,PM10的年平均濃度將平均上升0.262%,這說明第二產(chǎn)業(yè)的發(fā)展對PM10濃度帶來了重要的影響。第二產(chǎn)業(yè)主要分為工業(yè)和建筑業(yè),對于工業(yè)來說,其對煤炭、原油等化石燃料的燃燒是該行業(yè)產(chǎn)生PM10的主要原因。僅2013年,我國工業(yè)煤炭消費總量占煤炭消費總量的比重就已經(jīng)達(dá)到94.99%。因此,第二產(chǎn)業(yè)工業(yè)越發(fā)達(dá)的地區(qū),其PM10將可能越高。而建筑業(yè)由于對鋼筋、水泥的需求比較大,因此,也間接地增加了化石燃料的燃燒,另外,建筑施工過程引起的揚塵也是PM10的主要來源。從天然氣的供氣量(X10)來看,當(dāng)天然氣供氣總量每增加1%時,PM10的年平均濃度將平均上升0.296%,這結(jié)果可能與天然氣作為一種清潔能源相悖,但從目前我國天然氣的用途來看,其用途主要分為城市燃?xì)狻⒒ぴ稀⑻烊粴獍l(fā)電、工業(yè)燃料四大塊,并且化工原料和工業(yè)燃料的用氣量占絕大部分。2013年,化工原料和工業(yè)燃料的天然氣用氣量占天然氣用氣總量的41%。因此,天然氣供氣總量很大一部分還是用于工業(yè)發(fā)展,天然氣供氣總量大的地區(qū),其工業(yè)的比重也可能比較高,因此,工業(yè)污染對PM10的年平均濃度的影響可能也比較大。而森林覆蓋率(X16)與PM10的年平均濃度則呈負(fù)相關(guān)作用,當(dāng)森林覆蓋率每增加1%時,PM10的年平均濃度將平均下降0.308%,這主要是由于森林對PM10具有滯留、吸附和攔截功能,因此,森林覆蓋率越大,PM10的年平均濃度將越低。2.5.2 SO2、NO2的空間回歸分析與檢驗

表4 PM10的回歸分析與診斷性檢驗

同理,對SO2、NO2進(jìn)行空間回歸分析,得到下頁表5、表6。結(jié)果表明,SO2的空間回歸模型更適合空間滯后模型,其表達(dá)式如式(15)所示;而NO2則沒有通過拉格朗日乘數(shù)檢驗,其回歸模型使用普通最小二乘估計即可,將不顯著變量剔除后,得到修正后的模型結(jié)果如式(16)、下頁表7所示。從SO2模型中各變量的回歸系數(shù)和顯著性檢驗來看,在5%的顯著性水平下,各變量均通過了顯著性檢驗。在空間區(qū)域影響方面,空間回歸系數(shù)ρ對應(yīng)的P值為0.00055,說明我國各省份之間的SO2也存在空間溢出效應(yīng),并且當(dāng)空間滯后值每增加1%時,本地的SO2年平均濃度將平均增加0.491%。從社會發(fā)展影響來看,SO2年平均濃度與第二產(chǎn)業(yè)增加值占GDP的比重呈正相關(guān),當(dāng)比重每增加1%時,SO2年平均濃度將平均上升0.429%,這說明第二產(chǎn)業(yè)增加值所占地區(qū)經(jīng)濟(jì)的比重越大,其SO2年平均濃度將可能越高,這主要還是由于第二產(chǎn)業(yè)中的工業(yè)廢氣的排放巨大。2014年,全國廢氣中SO2排放量為1974.4萬噸,其中,工業(yè)SO2排放量為1740.4萬噸,所占比重達(dá)到88.15%。因此,治理工業(yè)廢棄項目對霧霾治理至關(guān)重要。從治理廢棄項目完成投資(X18)的回歸系數(shù)可以發(fā)現(xiàn),當(dāng)治理廢棄項目完成投資每增加1%時,SO2年平均濃度將平均下降0.270%。另外,森林也有降低SO2濃度的作用,當(dāng)森林覆蓋率每增加1%時,SO2年平均濃度將平均下降0.233%。這主要是由于植物對SO2有吸收作用,可對空氣進(jìn)行凈化。而對于NO2,X3(城鎮(zhèn)人口比重)和X21(民用汽車擁有量)與NO2年平均濃度呈正相關(guān),城鎮(zhèn)人口比重每增加1%,NO2年平均濃度將平均增加0.365%,民用汽車擁有量每增加1%,NO2年平均濃度將平均增加0.503%。從這兩個指標(biāo)可以看出,我國各省份大氣中的NO2濃度主要源自城市中汽車尾氣的排放,當(dāng)城市人口比重越大,所擁有的民用汽車數(shù)量將越多,其大氣中的NO2濃度將越高。此外,與PM10和SO2相似,提高森林覆蓋率將是降低NO2濃度一個有效的措施。

表5 SO2的回歸分析與診斷性檢驗

表6 NO2的回歸分析與診斷性檢驗

表7 NO2修正后的回歸模型

3 結(jié)論

本文通過空間計量分析和LASSO變量選擇方法,對中國31個省份造成霧霾天氣的三種主要污染物——PM10、SO2和NO2進(jìn)行研究分析。結(jié)果發(fā)現(xiàn):(1)華北地區(qū)及其周邊霧霾天氣最嚴(yán)重,呈“高高集聚”的空間分布。(2)第二產(chǎn)業(yè)是空氣中PM10和SO2的主要來源,且空間效應(yīng)顯著。第二產(chǎn)業(yè)中煤炭等化石燃料的燃燒和廢棄項目的排放是PM10和SO2的主要來源,通過治理工業(yè)廢棄項目完成投資的增加有助于SO2排放的減少。另外,各省份的PM10和SO2確實存在空間效應(yīng),采用空間滯后模型最合適。(3)汽車尾氣排放是NO2的主要來源。隨著城市化進(jìn)程的加快,人們的生活水平不斷提高,越來越多的人擁有自己的私人汽車,而隨著汽車數(shù)量的急劇增加,汽車的尾氣排放成為NO2的主要來源,并且其回歸模型采用最小二乘估計比較合適。

猜你喜歡
分析模型
一半模型
隱蔽失效適航要求符合性驗證分析
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動化發(fā)展趨勢分析
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
中西醫(yī)結(jié)合治療抑郁癥100例分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 亚洲精品无码抽插日韩| 欧美激情网址| av尤物免费在线观看| 成·人免费午夜无码视频在线观看| 国产在线观看精品| 欧美综合中文字幕久久| 伊人成人在线视频| 无码aaa视频| 亚洲AⅤ综合在线欧美一区| 国产欧美亚洲精品第3页在线| 亚洲色成人www在线观看| 一本大道视频精品人妻 | 综1合AV在线播放| 亚洲国产精品日韩欧美一区| 99视频有精品视频免费观看| 精品国产电影久久九九| 久久国产精品影院| 日韩麻豆小视频| 国产一二三区视频| 欧美日韩成人在线观看| 天天色综合4| 制服丝袜亚洲| www.91中文字幕| 亚洲中文字幕日产无码2021| 欧美在线伊人| 国产高清不卡| 99re在线免费视频| 亚洲国产综合精品一区| 五月天福利视频| 秋霞国产在线| 99视频有精品视频免费观看| 国产成熟女人性满足视频| 亚洲国产欧美目韩成人综合| 高清乱码精品福利在线视频| 国产成人精品视频一区视频二区| 在线视频亚洲欧美| 亚洲AV无码一区二区三区牲色| 污污网站在线观看| 国产女人综合久久精品视| 国产欧美日韩va另类在线播放 | 日韩高清中文字幕| 中文字幕 91| 波多野结衣第一页| 色精品视频| 国产成人区在线观看视频| 精品一区二区三区视频免费观看| 无码视频国产精品一区二区| 久久久91人妻无码精品蜜桃HD| 亚洲日本在线免费观看| 亚洲欧美日韩精品专区| 58av国产精品| 国产第二十一页| 国产地址二永久伊甸园| 日韩 欧美 国产 精品 综合| 国产精品第一区| 大陆精大陆国产国语精品1024| 国产乱子伦精品视频| 欧美午夜在线视频| 精品视频在线一区| 亚洲区一区| 亚洲高清在线播放| 亚洲第一成年免费网站| 国语少妇高潮| 精品亚洲麻豆1区2区3区| 91尤物国产尤物福利在线| 国产精品九九视频| 欧美精品在线看| 免费一级毛片在线观看| 国产不卡网| 91精品啪在线观看国产60岁| 国产麻豆精品在线观看| 成人伊人色一区二区三区| 亚洲三级电影在线播放| 亚洲日韩精品欧美中文字幕 | 中文字幕在线欧美| 欧美日韩激情在线| 日韩一区精品视频一区二区| 国产欧美日韩综合在线第一| 久热re国产手机在线观看| 亚洲无限乱码| 无码丝袜人妻| 久草热视频在线|