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

HBV水文模型在玉樹巴塘河流域洪水臨界雨量閾值研究中的應(yīng)用

2015-04-10 10:53:24劉義花魯延榮李昌玉
水土保持研究 2015年2期
關(guān)鍵詞:模型

劉義花, 魯延榮, 周 強(qiáng), 李昌玉

(1.青海省氣候中心, 西寧 810001; 2.青海省基礎(chǔ)地理信息中心, 西寧810000;3.青海省師范大學(xué) 生地學(xué)院, 西寧801600; 4.西寧市氣象局 西寧 810001)

?

HBV水文模型在玉樹巴塘河流域洪水臨界雨量閾值研究中的應(yīng)用

劉義花1, 魯延榮2, 周 強(qiáng)3, 李昌玉4

(1.青海省氣候中心, 西寧 810001; 2.青海省基礎(chǔ)地理信息中心, 西寧810000;3.青海省師范大學(xué) 生地學(xué)院, 西寧801600; 4.西寧市氣象局 西寧 810001)

區(qū)域氣象災(zāi)害的評(píng)估在防災(zāi)減災(zāi)中具有很重要的地位,它不僅是認(rèn)識(shí)災(zāi)情、進(jìn)行災(zāi)害區(qū)劃、實(shí)行災(zāi)害預(yù)測(cè)、制定防治對(duì)策、進(jìn)行損失評(píng)估、實(shí)施防治措施和進(jìn)行項(xiàng)目管理的基礎(chǔ),對(duì)政府的輔助決策都具有重要意義。基于玉樹縣社會(huì)經(jīng)濟(jì)統(tǒng)計(jì)資料、水文資料、巴塘河洪水災(zāi)情資料的基礎(chǔ)上,應(yīng)用HBV模型嘗試性的研究暴雨誘發(fā)的中小河流洪水臨界風(fēng)險(xiǎn)雨量閾值研究,結(jié)果表明:1)近11 a來新寨站平均流量23.2 m3/s,2001年、2003年、2005年汛期流量較大,2006—2011年流量明顯偏少;2)在率定期HBV模型對(duì)新寨站日徑流深模擬的確定性系數(shù)達(dá)0.678 2,Nash效率系數(shù)為0.604 4,驗(yàn)證期確定性系數(shù)超過了0.770,Nash效率系數(shù)為0.530 5;3)根據(jù)不同的基礎(chǔ)水位,有效劃分了24 h玉樹巴塘河流域洪水面雨量預(yù)警指標(biāo),為今后玉樹縣巴塘河流域提高災(zāi)害防御能力提供研究基礎(chǔ)。

巴塘河流域; HBV水文模型; 洪水

IPCC第四次評(píng)估報(bào)告指出[1],全球變暖導(dǎo)致極端事件頻發(fā),雖然極端事件是小概率事件,但是一旦發(fā)生將對(duì)自然和社會(huì)產(chǎn)生嚴(yán)重影響。在全球氣候變暖的大背景下,青海牧區(qū)表現(xiàn)出溫度升高、降水變率加大的區(qū)域響應(yīng),造成干旱、雪災(zāi)、暴雨洪澇等極端天氣氣候事件、氣象災(zāi)害加劇,青藏高原是氣候和生態(tài)環(huán)境變化的敏感區(qū)域和脆弱地帶,通過加強(qiáng)青藏高原氣象災(zāi)害特征和氣象災(zāi)害風(fēng)險(xiǎn)技術(shù)的研究,可提高對(duì)氣象災(zāi)害、災(zāi)害風(fēng)險(xiǎn)意識(shí),提高災(zāi)害預(yù)測(cè)預(yù)報(bào)水平,為政府和公眾提供準(zhǔn)確的預(yù)報(bào)和服務(wù),為制定農(nóng)牧業(yè)發(fā)展戰(zhàn)略,進(jìn)行正確的農(nóng)牧業(yè)決策,采取有效措施抗御自然災(zāi)害與保持農(nóng)牧業(yè)的高產(chǎn)和穩(wěn)產(chǎn)提供理論依據(jù)。根據(jù)1984—2007年青海省氣象災(zāi)害統(tǒng)計(jì),冰雹、干旱災(zāi)害年均造成農(nóng)業(yè)作物受災(zāi)面積400 km2以上,暴雨洪澇、雪災(zāi)、霜凍災(zāi)害所造成的農(nóng)業(yè)作物受災(zāi)面積也較大,分別為52.59,41.67,85.87 km2[2]。因此,通過各種技術(shù)和方法有效提高中小河流洪水、山洪地質(zhì)災(zāi)害風(fēng)險(xiǎn)評(píng)估和風(fēng)險(xiǎn)區(qū)劃能力,為青海省暴雨誘發(fā)的中小河流洪水、山洪地質(zhì)災(zāi)害氣象風(fēng)險(xiǎn)預(yù)警服務(wù)提供科技支撐。

國內(nèi)外針對(duì)臨界雨量閾值方面的研究主要有:何思為等[3]應(yīng)用四個(gè)水文模型對(duì)黑河流域上游進(jìn)行了研究對(duì)比,張建新等[4]應(yīng)用HBV模型開展對(duì)東北地區(qū)多冰雪地區(qū)洪水預(yù)報(bào)的研究,段生榮[5]以黃河流域大通河支流為例應(yīng)用統(tǒng)計(jì)方法研究臨界面雨量,張世才等[6]通過實(shí)際災(zāi)情研究祁連山區(qū)山洪臨界雨量并對(duì)該地區(qū)做了風(fēng)險(xiǎn)區(qū)劃。目前國內(nèi)外學(xué)者對(duì)青藏高原的生態(tài)環(huán)境變化做了大量研究,但是在氣象災(zāi)害風(fēng)險(xiǎn)評(píng)估方面的研究處于起步階段。筆者查閱大量文獻(xiàn)[7-11],發(fā)現(xiàn)針對(duì)青藏高原暴雨誘發(fā)的中小河流洪水臨界雨量研究較少,多采用統(tǒng)計(jì)方法,本研究嘗試性的采用HBV模型,通過對(duì)模型輸入?yún)^(qū)域水文資料、土地類型資料、氣象資料,從而直觀對(duì)比水文站觀測(cè)數(shù)據(jù)與模型模擬數(shù)據(jù)之間的差異,旨在為青海省玉樹巴塘河流域洪水預(yù)警提供具有客觀、科學(xué)性的臨界閾值,為今后巴塘河流域風(fēng)險(xiǎn)預(yù)警提供科學(xué)技術(shù)支撐。

1 研究區(qū)概況

巴塘河為通天河右岸一級(jí)支流[7],位于青海省玉樹藏族自治州玉樹縣境內(nèi),因流經(jīng)巴塘盆地而得名,巴塘河發(fā)源于格拉山北日阿如東塞以東4 km處,源流向北穿過峽谷進(jìn)入巴塘盆地,納支流扎巴曲后始稱巴塘河,后又進(jìn)入峽谷北上至玉樹縣城結(jié)古鎮(zhèn),最終匯入通天河。巴塘河全長(zhǎng)約92 km,流域面積1 793.9 km2,海拔2 778~5 672 m,徑流補(bǔ)給以降水為主,玉樹縣年平均氣溫3.8℃,年降水量為484.6 mm,流域內(nèi)土地類型以林地、草地、裸巖石礫地為主。從天氣學(xué)角度來看,巴塘河流域暴雨以中小尺度天氣系統(tǒng)為主,降水過程一般歷時(shí)短,強(qiáng)度大,再加上地形起伏度較大,因此一旦形成洪水,對(duì)農(nóng)牧業(yè)設(shè)施以及人類的生命財(cái)產(chǎn)安全帶來嚴(yán)重威脅。

2 研究方法簡(jiǎn)介

2.1 模型原理

HBV水文模型是瑞典國家水文氣象局開發(fā)研制,基于DEM劃分子流域的半分布式的概念性水文模型。多種情況下模型誤差小于20%。模型很簡(jiǎn)單,非常適用于大流域。模型不同版本已在全世界40多個(gè)位于不同氣候區(qū)的國家,如瑞典、津巴布韋、印度、哥倫比亞和中國等國家的洪水預(yù)報(bào)、水資源評(píng)估、營(yíng)養(yǎng)鹽負(fù)荷估算等領(lǐng)域得到廣泛應(yīng)用。HBV-D模型,由德國PIK研究所Krysanova.V博士改進(jìn)。HBV-D模型由氣候資料插值、積雪和融化、蒸散發(fā)估算、土壤濕度計(jì)算過程、產(chǎn)流過程、匯流過程等子模塊組成。該版本模型具有匯流時(shí)間模塊,分別模擬各子流域的徑流過程,后經(jīng)過河道匯流形成流域出口斷面的徑流過程,模型應(yīng)用相對(duì)簡(jiǎn)便,輸入數(shù)據(jù)主要是研究區(qū)DEM、日均氣溫、降雨、土地利用、土壤最大含水量和河流匯流時(shí)間等參數(shù),通過模型中實(shí)際徑流量與模擬值得比較,從而有效劃分了不同水位臨界高度下的洪水臨界雨量值,以此達(dá)到巴塘河流域洪水預(yù)警的目的。

2.2 風(fēng)險(xiǎn)雨量計(jì)算方法

流域內(nèi)發(fā)生不同等級(jí)的洪水,洪水匯集到河流中達(dá)到一定水位高度所對(duì)應(yīng)的某時(shí)段的降水量即為該淹沒等級(jí)的風(fēng)險(xiǎn)雨量,本研究中把水文站達(dá)到警戒水位、保證水位、堤防高度時(shí)的雨量分別做為三個(gè)等級(jí)的風(fēng)險(xiǎn)雨量。通過HBV模型模擬研究時(shí)間段徑流深度以及觀測(cè)流量和水位、水位和研究區(qū)面雨量的關(guān)系,劃分出研究區(qū)不同基礎(chǔ)水位下臨界面雨量的閾值。

2.3 水文要素的變化特征

新寨水文站是巴塘河的觀測(cè)水文站,現(xiàn)已收集到2001—2011年逐日水位和流量資料,因2007年以前的水位基面跟后期的不一致,因此,在建模的時(shí)候涉及基礎(chǔ)水位跟流量建立關(guān)系的時(shí)候,只考慮2007年水位和流量的變化關(guān)系。從新寨水文站流量圖上來看(圖1),近11 a來新寨站平均流量23.2 m3/s,每年汛期6—9月流量達(dá)到峰值,2001年、2003年、2005年汛期流量較大,直至2006—2011年流量明顯偏少,從水文站水位的變化來看,2011年汛期,有2次洪水過程達(dá)到警戒水位,有1次洪水過程達(dá)到保證水位,玉樹縣是半農(nóng)半牧的地區(qū),水文站的流量明顯變化與水文站的調(diào)蓄作用比較密切。

3 HBV模型的應(yīng)用

3.1 徑流深度模擬結(jié)果

利用水系、DEM等地理信息,采用GIS和水文分析技術(shù)提取了巴塘河流域的范圍、流域中心點(diǎn),基于R雨量插值軟件和流域內(nèi)氣象站觀測(cè)數(shù)據(jù),得出2000—2011年流域面雨量的逐日變化序列,因在每一次暴雨過程中,降雨徑流關(guān)系極為復(fù)雜,因此模型通過土壤持水力以及土壤類型、匯流到子流域出口的時(shí)間來共同模擬每一次洪水過程。在模型運(yùn)行之前,需要修改研究流域的面積,模擬的時(shí)間段,模型中涉及31個(gè)參數(shù),因此對(duì)每個(gè)參數(shù)敏感性分析具有重要意義,模型中相關(guān)參數(shù)每調(diào)整一次,模擬結(jié)果會(huì)發(fā)生明顯變化,因此用確定性系數(shù)和Nash效率系數(shù)驗(yàn)證模型的穩(wěn)定性和可靠性,該方法用來解釋模型的誤差,最終的目地是為了確定性系數(shù)和Nash效率系數(shù)達(dá)到0.5~0.9,以此來驗(yàn)證模型擬合的效果,如果模型中調(diào)參不合適,那么Nash效率系數(shù)是負(fù)值,模擬的徑流深度和實(shí)際觀測(cè)值的擬合效果偏差較大,基于此,對(duì)其中5個(gè)參數(shù)進(jìn)行了多次調(diào)整和敏感性分析(表1),將這些參數(shù)的不同組合輸入到HBV水文模型,率定后的HBV模型對(duì)新寨站日徑流深模擬的確定性系數(shù)達(dá)0.6782,Nash效率系數(shù)為0.6044(圖2a),尤其是2002年以來的模擬值與實(shí)際觀測(cè)值的擬合程度較好,能夠準(zhǔn)確捕捉每次洪水過程,模型模擬的結(jié)果與實(shí)況較一致,能夠很好地模擬出巴塘河流域的日徑流過程。為進(jìn)一步檢驗(yàn)HBV 模型效果,使用2007—2011年逐日氣象、水文資料對(duì)巴塘河流域的預(yù)報(bào)效果進(jìn)行了檢驗(yàn)(圖2b),可以看出經(jīng)過率定后的HBV 模型在巴塘河流域具有很強(qiáng)的適用性,對(duì)新寨站逐日徑流深模擬的確定性系數(shù)超過了0.770,Nash效率系數(shù)為0.530 5,模擬出的水文過程線與實(shí)際基本吻合,很好的預(yù)報(bào)出了洪水對(duì)降水的響應(yīng)過程,尤其是2007年、2010年、2011年以來每年的日徑流深度的擬合效果接近實(shí)際,以此為依據(jù),來確定洪水發(fā)生時(shí)達(dá)到不同高度的臨界面雨量是具有意義,并為科學(xué)、客觀預(yù)測(cè)強(qiáng)降水對(duì)當(dāng)?shù)卦斐珊蔚蕊L(fēng)險(xiǎn)提供科學(xué)技術(shù)支撐。

圖1 2007-2011年新寨水文站逐日水位和流量的變化

圖2 2000-2011年逐日模擬徑流深度與觀測(cè)徑流深度的對(duì)比

圖2中的nash系數(shù)和確定性系數(shù)是根據(jù)公式(1)和公式(2)計(jì)算得出的,如果nash系數(shù)和確定性系數(shù)不再0.5~0.9范圍內(nèi),則說明模擬值和觀測(cè)值擬合效果不好,模擬結(jié)果沒有實(shí)際意義。

(1)

(2)

3.2 水位和流量的關(guān)系

通過水文站逐日觀測(cè)流量和水位,建立2007—2011年洪水流量和水位之間的關(guān)系,流量和水位相關(guān)性特別高,通過0.001的顯著性檢驗(yàn)。由圖3可知,自2007年以來,流量和水位呈現(xiàn)持續(xù)上升的態(tài)勢(shì),1971—2000年青南牧區(qū)年平均降水量480.1 mm,2001—2010年年平均降水量499.7 mm,說明青南牧區(qū)降水明顯增加,降水量增加趨勢(shì)與巴塘河流量增加的趨勢(shì)一致。通過HBV模型模擬結(jié)果以及流量和水位的換算關(guān)系,可以得出逐日模擬流量、模擬水位以及流域面雨量,為流域臨界面雨量的閾值劃分提供基礎(chǔ)。

表1 參數(shù)的意義以及敏感性

4 臨界雨量閾值的劃分

根據(jù)玉樹新寨站警戒水位、保證水位、堤壩高度,根據(jù)洪水達(dá)到警戒、保證或漫過堤壩水位的不同條件來判定暴雨誘發(fā)的洪水對(duì)整個(gè)巴塘河流域造成不同的風(fēng)險(xiǎn)等級(jí)。因高原獨(dú)特氣候條件,通過分析1961—2011年歷史洪水災(zāi)情以及水位、流量的變化情況,發(fā)現(xiàn)1 h,3 h,6 h,12 h的暴雨的強(qiáng)度不是特別大,因此劃分24 h不同基礎(chǔ)水位下的山洪致災(zāi)臨界雨量是符合客觀實(shí)際,新寨站警戒水位、保證水位、堤防高度分別為3 655.77,3 655.83,3 656.91 m,研究達(dá)到不同水位高度時(shí)的臨界面雨量(圖4),劃分出臨界面雨量對(duì)每次洪水過程有較好的響應(yīng),氣象部分根據(jù)面雨量預(yù)報(bào)值判定降水量是否造成洪水,從而達(dá)到洪水風(fēng)險(xiǎn)預(yù)警的目的。

圖3 2007-2011年逐日流量與水位擬合

圖4 不同基礎(chǔ)水位下臨界面雨量值

5 結(jié)論與討論

本研究基于巴塘河流域高程圖、土地類型、氣溫、降水以及土壤持水量等要素模擬巴塘河流域模擬徑流深度,率定期的HBV 模型對(duì)新寨站日徑流深模擬的確定性系數(shù)達(dá)0.678 2,驗(yàn)證期確定性系數(shù)超過了0.770 0,模擬出的水文過程線與實(shí)際基本吻合,很好的預(yù)報(bào)出了洪水對(duì)降水的響應(yīng)過程。通過觀測(cè)流量和水位擬合關(guān)系以及流量與面雨量的關(guān)系,以洪水達(dá)到警戒、保證或漫過堤壩水位的不同條件來判定暴雨誘發(fā)的洪水對(duì)整個(gè)巴塘河流域造成不同的風(fēng)險(xiǎn)等級(jí)。通過半分布式水文模型在青藏高原暴雨誘發(fā)的中小河流洪水臨界面雨量研究中的應(yīng)用,氣象部門通過該流域臨界雨量閾值的范圍,如果該流域面雨量達(dá)到一定的量級(jí),可有效、快速的發(fā)布預(yù)警,為風(fēng)險(xiǎn)預(yù)警提供科學(xué)依據(jù)。由于高原本身的地域特色,氣象數(shù)據(jù)、水文數(shù)據(jù)、災(zāi)情數(shù)據(jù)的不完備,短期內(nèi)無法驗(yàn)證劃分閾值的有效性,隨著今后汛期強(qiáng)降水過程,逐步驗(yàn)證模型劃分臨界面雨量閾值的有效性并不斷完善。

[1] Solomon S, Qin D, Manning M, et al. Climate Change 2007:The physical Science Basis:Contribution of working group I to the Fourth Assessment Report of the intergovernmental Panel on climate change. IPCC Fourth assessment Report [M]. Cambridge University Press,2007.

[2] 劉義花,李林,蘇建軍,等.青海省春小麥干旱災(zāi)害風(fēng)險(xiǎn)評(píng)估與區(qū)劃[J].冰川動(dòng)土,2012,34(6):1416-1423.

[3] 何思為,南卓銅,王書功,等.四個(gè)概念水文模型在黑河流域上游的應(yīng)用與對(duì)比分析[J].水文,2012,32(3):13-19.

[4] 張建新,趙孟芹,章樹安,等.HBV模型在東北多冰雪地區(qū)的應(yīng)用研究[J].水文,2007,27(4):31-33.

[5] 段生榮.典型小流域山洪災(zāi)害臨界雨量計(jì)算分析應(yīng)用研究[J],中國農(nóng)村水利水電,2008(8):63-68.

[6] 張世才,褚建華,張同澤,等.祁連山區(qū)山洪災(zāi)害臨界雨量計(jì)算分析和風(fēng)險(xiǎn)區(qū)劃[J].水土保持學(xué)報(bào),2007,21(5):196-200.

[7] 李燕萍.玉樹巴塘河初步水文分析與計(jì)算[J].青海科技,2010(4):47-50.

[8] 趙彥增,張建新,章樹安,等,HBV模型在淮河官寨流域的應(yīng)用研究[J].水文,2007,27(2):57-59.

[9] 祁元,劉勇,楊正華,等.基于GIS蘭州滑波與泥石流災(zāi)害危險(xiǎn)性分析[J].冰川凍土,2012,34(1):96-104.

[10] 文明章,林昕,游立軍,等.山洪災(zāi)害風(fēng)險(xiǎn)雨量評(píng)估方法研究[J].氣象,2013,39(10):1325-1330.

[11] 陳曉弟,羅京義,謝仁波,等.銅仁錦江河流域面雨量計(jì)算方法探討[J].貴州氣象,2010(S):134-137.

Application of HBV Model to the Study on Risk Precipitation in Different Grades in Batang River Region

LIU Yihua1, LU Yanrong2, ZHOU Qiang3, LI Changyu4

(1.ClimateCenterofQinghai,Xi′ning810001,China; 2.QinghaiProvinceBasicGeographicInformationCenter,Xi′ning810000,China; 3.CollegeofLifeScience,QinghaiNormalUniversity,Xi′ning801600,China; 4.Xi′ningMeteorologyBureau,Xi′ning810001,China)

The trend of extreme climate events has been analyzed by many researchers. Once it occurs, it maybe cause serious influence, so it becomes more and more important to research meteorology disaster, which can not only help us know disaster, conduct disaster forecast and loss assessment, but also help government to solve the problems. The module of runoff was simulated using a hydrological model(HBV) based on meteorological data, land use, hydrological data and Dem of Batang river region. The results of simulated runoff depth were verified. The results showed that: 1) the average runoff is 23.2 m3/s in recent 11 years in Xinzhai station, the runoff volumes in flood season in 2001,2003 and 2005 were greater, but the runoff volumes became less from 2006 to 2011; 2) the coefficients of determination and NASH were 0.678 2, 0.604 4, respectively, in the calibrated period, and were 0.770 and 0.530 5, respectively, in the verification period; 3) the warning critical rainfalls triggering runoff at different water levels were calculated using HBV model in order to reduce losses resulting from the flooding disasters.

Batang river region; HBV model; mudflow

2014-06-09

2014-07-10

國家自然科學(xué)基金項(xiàng)目“湟水流域農(nóng)牧民對(duì)氣候變化的適應(yīng)行為與感知基礎(chǔ)研究”(41261010)

劉義花(1979—),女,青海平安人,碩士,主要從事氣候評(píng)價(jià)方面工作。E-mail:yihualiu12@126.com

P338

1005-3409(2015)02-0224-05

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产在线精品美女观看| 国产性爱网站| 亚洲人成网站观看在线观看| 久久综合一个色综合网| 日韩天堂在线观看| 国产女人18水真多毛片18精品| 成人午夜视频网站| 四虎影院国产| 亚洲乱码在线视频| 伊人色天堂| 无码在线激情片| 国产成年无码AⅤ片在线| 91精品国产丝袜| 无码一区18禁| 日本一本正道综合久久dvd| 国产视频只有无码精品| 国内精品91| 国产欧美亚洲精品第3页在线| 国产成人高清精品免费软件| 日韩久久精品无码aV| 美女啪啪无遮挡| 亚洲熟女偷拍| 国产麻豆精品久久一二三| a网站在线观看| 动漫精品啪啪一区二区三区| 色哟哟国产精品一区二区| 久久这里只有精品国产99| 国产成人精品一区二区免费看京| 国产91久久久久久| 国产网站一区二区三区| 成人年鲁鲁在线观看视频| 亚洲成a∧人片在线观看无码| 中文字幕欧美成人免费| a在线亚洲男人的天堂试看| 亚洲欧美另类日本| 无码专区国产精品第一页| 日韩欧美在线观看| 国产精品夜夜嗨视频免费视频| 伊人久久大香线蕉成人综合网| 无码视频国产精品一区二区| 青青草原国产| 日韩色图区| 久久综合丝袜长腿丝袜| 国产黄色爱视频| 色国产视频| 欧美一区二区三区香蕉视| 久久九九热视频| 日本不卡在线播放| 国产99欧美精品久久精品久久 | 日本五区在线不卡精品| 国产波多野结衣中文在线播放 | 国产精品开放后亚洲| 日a本亚洲中文在线观看| 毛片久久网站小视频| 无码在线激情片| 日本AⅤ精品一区二区三区日| 欧美不卡二区| 久久99国产综合精品1| 欧美成人区| 午夜啪啪网| 久久五月天综合| 久久中文字幕2021精品| 亚洲高清无在码在线无弹窗| 国产一区免费在线观看| 国产精品无码AV中文| 亚洲av日韩av制服丝袜| 欧美.成人.综合在线| 亚洲一级色| 日韩精品专区免费无码aⅴ| 99精品伊人久久久大香线蕉| 91 九色视频丝袜| 中文字幕亚洲综久久2021| 草草影院国产第一页| 97国产精品视频自在拍| 亚洲三级色| 乱人伦视频中文字幕在线| 国产精品人莉莉成在线播放| 亚洲三级成人| 精品黑人一区二区三区| 在线免费不卡视频| 成人一级免费视频| 色综合中文综合网|