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

水文時間序列高階自相關性識別及其對水文分析的影響研究

2019-02-28 06:31:32,,2,,
人民珠江 2019年2期
關鍵詞:趨勢

,,2,,

(1. 長安大學 環境科學與工程學院,陜西西安710054;2. 長安大學 旱區地下水文與生態效應教育部重點實驗室,陜西西安710054;3. 陜西省江河水庫管理局,陜西西安710018)

水文時間序列是指某種水文特征值隨時間而變的一系列觀測值,在社會水文學的視點下,它也可以被看作是由氣候條件、下墊面條件和人類活動等因素綜合作用結果的表征[1]。一般水文時間序列可分解為確定性和隨機性兩種成分,確定性成分主要表現為水文現象的趨勢變化和周期變化等,其主要受人類活動、氣候因素和下墊面因素的宏觀尺度驅動,使水文時間序列產生緩慢的漸變或劇烈的突變,亦或有規律的波動;而隨機性成分則表現為水文現象的相依性和純隨機變化,其主要受氣候、下墊面等各種隨機因素影響,并驅動水文時間序列的隨機性變化和不確定性變化[2]。

近年來,受全球氣候變化及高強度人類活動(如水庫蓄水、淤地壩攔水、引水灌溉及城鄉建設等)的影響,流域的水循環過程和水文物理形成機制都發生了顯著的改變[3-4]。很多地區的實測年徑流、降水、泥沙等水文時間序列都表現出相對較大的自相關性[5-6],使得一些檢驗結果的準確性受到一定影響。目前,關于水文序列自相關性的研究工作主要集中在序列的一階自相關性上。如在國內外相關領域廣泛應用的由Von Storch和NaVarra提出的預置白[7](Pre-whitening,PW)方法和Yue等提出的去趨勢預置白[8](Trend-free pre-whitening,TFPW)方法等,其主要作用都是消除或降低序列中的一階自相關性,以滿足數據獨立性的要求。然而,RazaVi和Vogel[9]的研究表明,在對樹木年輪序列進行預置白(Pre-whitening,PW)處理后會使時間序列中包含的重要信息丟失,同時會導致對氣候和水文極端條件的低估,并建議在使用預置白方法時需要考慮其在長期自相關性上的陷阱。但對于序列高階自相關性的剔除方法或影響,這些文獻均未有過多探討。目前,水文領域有關高階自相關問題的論述,并不多見。張洪波等[10]提出基于EMD的改進RBF方法來控制高頻分量的預測誤差,進而實現提高徑流預測精度的方法。其中作者根據各分量的自相關分析結果,篩選存在顯著自相關的時滯,完成選階后構建神經網絡的輸入向量,盡管考慮了各顯著階數序列,但未說明其對徑流預測效果產生的影響。Sen[11]提出的一種新的用于統計檢驗的自相關處理方法,即過白化處理法(over-whitening,OW),其通過在時間序列的標準化序列上添加一個白噪聲,可實現一階和高階自相關性的綜合去除,為高階自相關移除提供了手段。此外,Julian[12]對德國1978—2009年17個觀測站的降水序列中18O、2H穩定同位素進行趨勢變化研究,結果表明在建模過程中考慮高階自相關性的影響時,有3個站的序列表現出顯著的趨勢變化,而廣泛應用的僅考慮一階自相關性的趨勢分析方法未能充分考慮到穩定同位素序列中高階自相關性的存在。再者,張洪波等[13]在傳統過白化處理法的基礎上,提出了改進措施,采用分段過白化的處理手段實現了趨勢低損情況下的高階顯著自相關剔除。綜上所述,在對水文序列進行分析時,僅僅考慮序列的一階自相關性是不夠充分的。與此同時,盡管水文序列的高階自相關議題在一些文獻中出現過,但其并未引起廣泛的關注或普遍的共識。因此,非常有必要對水文序列中高階自相關性做以深入的解析并探討其對現有技術手段的影響。

為此,本文分別選取了陜西境內漢江流域的2個水文站、渭河流域的6個水文站以及無定河流域的11個水文站的實測年徑流序列為研究對象,通過分段去趨勢方法,剔除序列中的趨勢成分,再對去趨勢序列進行自相關性分析,以確定序列高階自相關性在不同流域的空間分布。同時,文章進一步以BP神經網絡為基礎,建立降水-徑流預測模型,旨在探討加入高階降雨序列對徑流預測效果的影響,進而提出水文時間序列分析中高階自相關性的處理方法。

1 研究區概況與數據來源

1.1 研究區概況

1.1.1漢江流域

漢江是陜南地區最大的河流,也是長江較大的一級支流,發源于陜西省寧強縣潘家山,由西向東流經漢中、安康兩市進入湖北省境。漢江在陜西境內河長654 km,集水面積4.86萬km2,多年平均徑流量247億m3。徑流量地區分布不均勻,總的趨勢是南岸多于北岸,即大巴山多,秦嶺少。陜南地區降水充沛,時空分布不均,年均降水量1 050 mm,秦嶺南坡年均降水量在600~800 mm之間,巴山北坡米倉山區年降水量最大,可達1 400 mm以上。流域內河川徑流的補給主要以降雨形成的地表徑流為主,地下水補給為次[16]。

漢江沿岸的地貌特征以丘陵、盆地為主,東西長約400 km,南北寬約3~60 km,漢江流貫中部,形成以漢中盆地和石泉-安康盆地為主的多個盆地區[17]。

1.1.2渭河流域

渭河是黃河的一級支流,起源于甘肅省渭源縣,流經甘肅、寧夏、陜西,于潼關匯入黃河,流域總面積約13.43萬km2,全長818 km,陜西境內河流總長約502 km,流域面積6.71萬km2。渭河流域年平均氣溫12 ℃~14 ℃,年平均降雨量為400~800 mm,降雨主要集中在6—9月,降雨量占全年的60%左右[13-14]。

渭河流域地形總特點是西北高東南低,陡坡自西向東逐漸變緩,河谷自上游至下游逐漸變寬。渭河流域北岸支流多發源于黃土丘陵和黃土高原,相對源遠流長,比降較小,含沙量大;南岸支流均發源于秦嶺山區,源短流急,谷狹坡陡,徑流較豐,含沙量小[15]。

1.1.3無定河流域

無定河是榆林地區最大的河流,它發源于定邊縣白于山北麓,流經米脂、綏德到清澗縣川口以南20 km處注入黃河。全長491.0 km,流域面積30 260 km2,陜西境內河長442.8 km,流域面積21 049.3 km2,流域降水量從東南向西北遞減,多年平均降雨量為409 mm[10]。

按地形地貌及水土流失特點,流域可分為風積沙丘沙蓋區、黃土斜梁高山區和黃土丘陵溝壑梁峁區三大區域。流域西北部為風積沙丘沙蓋區,約占全面積的54.3%,氣候干燥,地面起伏平緩,水網稀少;西南部為黃土斜梁高山區,約占全面積的11.4%,地面呈梁峁狀丘陵;東南部為黃土丘陵溝壑梁峁區,約占全面積的34.3%,地表支離破碎,梁峁起伏,溝壑縱橫。3個流域的自然地理條件對比見表1。

注:表1相關統計數據出自于陜西省2016年水資源公報

1.2 數據來源

本研究采用的數據主要有漢江流域石泉、安康水文站,渭河流域華縣、臨潼、咸陽、魏家堡、林家村、北道水文站以及無定河流域韓家峁、橫山、趙石窯、白家川、殿市、馬湖峪、青陽岔、李家河、曹坪、丁家溝、綏德水文站等共19個水文站點的實測年徑流及控制流域的面降雨數據。其中,年徑流數據來源于黃委水文資料匯編,降水數據來源中國氣象資料共享網,站點數據經加權平均后獲得面降雨量數據。各流域站點分布見圖1。

圖1 各流域站點分布

2 研究方法

2.1 分段去趨勢

受水文序列非一致性的影響,其本身表現出時域上的差異性特征,而非單調性的趨勢變化。因此,在進行自相關性分析前,需要對水文序列進行去趨勢處理。本文使用了一種基于回歸變異的分段去趨勢方法來去除徑流序列的趨勢成分。其具體步驟如下:①采用Tomé和Miranda[18]提出的分段線性擬合技術,對水文時間序列進行回歸變異識別,得到該序列的多個回歸變異點;②擬合各分段回歸系數與回歸方程,分段去除趨勢,得到去趨勢序列,即平穩序列;③對去趨勢序列進行自相關性分析。

2.2 BP神經網絡

BP(back propagation)神經網絡是1986年由Rumelhart和McClelland為首的科學家提出的一種按照誤差逆向傳播算法訓練的多層前饋神經網絡,是目前應用最廣泛的神經網絡模型之一。BP網絡能學習和存貯大量的輸入-輸出模式映射關系,而無需事前揭示描述這種映射關系的數學方程。它的學習規則是使用最速下降法,通過反向傳播來不斷調整網絡的權值和閾值,使網絡的誤差平方和最小[19-20]。

本文中,BP神經網絡主要基于MATLAB軟件中的NEWFF函數開展,構建多輸入單輸出的神經網絡模型,其中訓練誤差GOAL設置為0.000 1;其他參數均采用默認值。經過反復迭代后,可確定各輸入變量的BP神經網絡模型,得到預測結果。

3 計算結果

3.1 高階自相關識別

3.1.1漢江流域

通過分段去趨勢方法對漢江流域的石泉和安康水文站的實測年徑流序列進行處理,得到其去趨勢序列,并采用偏自相關(PACF)函數對其進行自相關性檢驗。結果表明石泉、安康兩站點的去趨勢年徑流序列的偏自相關系數值均未超出臨界值,即不存在顯著的高階自相關性,見圖2。

3.1.2渭河流域

利用分段去趨勢方法對渭河干流6個水文站的年徑流序列進行處理,得到其去趨勢序列,并采用偏自相關(PACF)函數對其進行自相關性檢驗。結果與漢江流域類似,華縣、臨潼、咸陽、魏家堡、林家村、北道6個站點的去趨勢年徑流序列的偏自相關系數值均未超出臨界值,即均不存在高階自相關性,見圖3。

a) 石泉

b) 安康圖2 去趨勢年徑流序列PACF檢驗

a) 華縣

b) 臨潼

c) 咸陽

d) 魏家堡

e) 林家村

f) 北道圖3 去趨勢年徑流序列PACF檢驗

3.1.3無定河流域

將分段去趨勢方法應用于無定河流域11個水文站的年徑流序列,得到去趨勢序列,并采用偏自相關(PACF)函數對其進行自相關性檢驗。結果表明韓家峁站存在6階顯著高階自相關性;橫山站存在7階顯著高階自相關性,6階偏自相關系數值未超出臨界值,但數值較大,在后續計算中應考慮6階的影響;趙石窯站存在11階顯著高階自相關性;白家川站存在2階和5階顯著高階自相關性;殿市存在一階自相關;馬湖峪、青陽岔、李家河、曹坪、丁家溝、綏德此6個站點偏自相關系數值均未超出臨界值,即均不存在顯著自相關性,具體結果見圖4。

a) 韓家峁

b) 橫山

c) 趙石窯

d) 白家川

e) 殿市

f) 馬湖峪圖4 去趨勢年徑流序列的PACF檢驗

g) 青陽岔

h) 李家河

i) 曹坪

3.2 高階自相關的物理解析

從上節的高階自相關識別結果來看,不同的流域檢驗的結果差異較大。漢江流域和渭河流域研究站點均未檢測到高階自相關,而位于陜北的無定河流域則在多站點出現了高階自相關成分。一般高階自相關性主要受水循環因子的周期變化以及地下水入滲到排泄的時差的影響。由于并不是3個流域的水文序列都表現為高階自相關,因此,本節主要從地下水文循環的物理機制角度,分析產生這種空間差異的原因。

j) 丁家溝

k) 綏德續圖4 去趨勢年徑流序列的PACF檢驗

3.2.1漢江流域

在漢江流域選取的石泉和安康2個站點均位于安康市,區域地貌單元為安康盆地漢階地地層,其中淺層含水巖組的含水層巖性為粗砂、礫石、卵石等,深度在18~33 m以下,含水巖組受大氣降水及區域地下水補給,承壓性一般,水量較豐富;淺層承壓水含水巖組的含水層巖性為粗砂礫石、卵礫石、漂卵石等,深度在70~80 m以下,含水巖組受地下水補給,承壓性好,水量較豐富[21]。同時該區域雨量充沛,地表水豐富,為地下水形成提供了良好的補給條件,但地形破碎,河溝縱橫,儲水條件差,不利于地下水貯存[27]。因此,地下水補給地表徑流量的部分較少,且由于地處山區補給過程所需的時間跨度也較短,故徑流序列中產生高階自相關性的物理機制并不顯著存在,所以未能在序列中檢測到高階自相關性。

3.2.2渭河流域

渭河流域干流上、中、下游河水與地下水轉化關系有所不同。在渭河干流,上游地區地下水補給河水;中游地區河水一側補給地下水、另一側地下水補給河水,豐、枯季節地下水與河水的轉化關系也時常發生變化;下游地區如渭南至潼關段,河床淤積抬升,一般河水補給地下水[15]。但由于渭河徑流量中地下水補給量有限,加之渭河流域地處關中平原區,受當地水文地質條件影響,流域地下水循環較快,從入滲到排泄至河流的時滯較短,因此,其反映到年際尺度的水文序列變化上并不顯著。同時,近50 a來渭河流域降水量呈顯著減少趨勢[26],導致以降水為地下水主要補給源的渭河流域的地下水蓄水量補充有所削減,加之區內大量開采地下水,導致渭河基流量顯著減少,河川徑流中原有的地下水循環所蘊含的周期規律和相關關系在愈發弱化。由此可見,多因素的綜合影響導致了渭河流域各水文站的年徑流序列中未檢測到高階自相關性。

3.2.3無定河流域

無定河流域與漢江、渭河流域有較大不同,其河川徑流中30%~80%來自地下水補給[29],因此,地下水循環中的典型規律會較為顯著地體現在河川徑流過程中。

從地貌類型來看,無定河流域風沙區的韓家峁、橫山、趙石窯序列和黃土丘陵溝壑梁峁區的白家川序列檢驗出高階自相關性;其余站序列均未檢驗出高階自相關性。

無定河流域風沙區地勢平坦,林草覆蓋較好,土質為顆粒較大的沙土,地面滲漏強烈,下滲能力強[22],地下水埋深較淺,降水可補給到淺層地下水,而后再經由地下水徑流途徑,轉化為基流匯入河道。風沙區地面滲漏強烈,地下水補給徑流所占比重較大,一般達80%~90%以上[29],因此地下水循環特征可得到較好的表征。而黃土區植被條件差,地形起伏大、溝谷深切,黃土區地面下3.5 m處的恒濕層阻斷了地表水與地下水之間的水力聯系[25],使得地下水難以得到降水的有效補給,滲入土壤的水分通常無法到達飽和含水層,而以壤中流的方式補給到河流,因此在黃土區地下水循環路徑一般較短,短時滯在年際尺度的河川徑流變化中無法得到有效體現。黃土區地下水補給只占年徑流的30%左右[28-29]。白家川水文站是無定河流域的出口水文站,年徑流量遠大于其余水文站,且基徑比達71%[27];由于黃土區地下水補給量相對較小,因此其更趨向于表現上游主導地下水過程的相關或耦合特征。因此,上述因素可較為清楚地解釋風沙區年徑流序列呈現出高階自相關性,而黃土區的年徑流序列不存在高階自相關性的現象。

除此之外,筆者在研究中發現,出現高階自相關的站點所控制的流域都存在與其他流域存在側向交換的情況,即在某一含水層與其他流域存在地下水連通或補給關系。因此筆者認為無定河流域存在水文序列的高階自相關性可能與流域地下水分水嶺不閉合存在一定的關系。但限于目前無相關成果,無法佐證。

3.3 驗證分析

上文對韓家峁、橫山、趙石窯、白家川4個站點顯著高階自相關性的物理成因進行了簡單的定性分析,物理機制定量論證受觀測資料限制目前尚無法開展。本文擬以BP神經網絡為基礎構建降水-徑流模型,通過分析加入高階降雨序列對徑流預測效果的影響,從側面論證高階自相關的客觀存在性以及影響。

以橫山站實測年徑流和降雨序列(1957—2010年)為例,檢測結果顯示去趨勢年徑流序列存在7階顯著高階自相關性,同時6階偏自相關系數較大,不可忽視,故選取當年降雨序列P0、1—7階降雨序列(表示為P1…P7),當年去趨勢年徑流序列R0、進行BP神經網絡計算分析,模型的構建過程見圖5(P7+P0表示將兩組序列進行組合,其余類似)。

橫山站設置序列的訓練期為1964—1998年,輸入變量有5組,輸出變量有一組,將訓練期的輸入和輸出變量分別代入BP神經網絡,構建多輸入單輸出的BP神經網絡模型并得到5組訓練期預測值;設置序列的驗證期為1999—2010年,以構建的BP神經網絡模型為基礎,將5組輸入變量代入其中,得到5組驗證期預測值,并將其與5組訓練期預測值進行組合,得到的5組預測序列分別記為R1—R5,其與當年去趨勢年徑流序列R0的對比,見圖6。

圖5 BP神經網絡模型的構建過程(以橫山站為例)

由圖6和表2可以看出,當輸入變量依次為P0、P7+P0、P6+P7+P0、P1+P6+P7+P0即加入了不同階數的降雨序列時,預測序列R1—R4的預測效果逐漸變好,與R0的擬合度R2由0.336逐漸增加到0.769,表明當把不同階數的降雨序列加入當年降水序列中進行預測時能有效提高預測的效果,從側面驗證了徑流序列高階自相關的存在。但當加入1—6階的全部降雨序列時,預測序列R5的預測效果反而有所下降,與R0的擬合度降低至0.375。由此可見,在加入不同階數計算時應該首要考慮的是顯著階數、偏自相關系數較大值對應的階數,而把所有階數均加入計算時則有可能因為各階序列的重復疊加,使預測出現疊加誤差,進而導致預測效果有所下降。

a)對比R1

b) 對比R2

c) 對比R3

d) 對比R4

表2 橫山站5組預測序列與R0的擬合度

上文中所述的5組預測序列和R0均是不包含趨勢成分的殘余序列,并不存在實際的物理意義,遂將分段去趨勢步驟中所剔除的趨勢成分加回到預測序列與R0中,即得到重組序列RR1—RR5和實測年徑流序列R。此時重組序列則對應于實測年徑流序列的預測序列,重組序列與實測年徑流序列的相關關系,見圖7。與表2比較可知,將趨勢成分加回后,5組重組序列與實測年徑流序列的擬合度均有較大程度的提高。

a) RR1

b) RR2

c) RR3

d) RR4

e) RR5圖7 橫山站重組序列與實測年徑流序列的相關關系

將同樣的方法應用于韓家峁、白家川和趙石窯站,在構建BP神經網絡模型時所涉及的變量見表3。 計算得到3個站的預測序列與當年去趨勢年徑流序列、重組序列與實測年徑流序列的擬合度見表4。

由表4可知,當把顯著階數的降雨序列納入預報因子并實施預測時,可有效提高相關的徑流序列擬合度的大小;但當加入過多階降雨序列時,擬合度反而有所下降。基于以上分析可知,徑流序列中的高階自相關性是序列蘊含的規律性成分,在預測模型中考慮該特征,可有效提高徑流預測的效果。由此,在僅有短期實測徑流資料情況下,需要展延年徑流資料時,可利用序列的顯著低階與高階自相關性,優化輸入變量,構建多輸入單輸出的BP神經網絡模型,可在一定程度上有效提高預測的效果,并不失為一種提高徑流預測效果的新思路。

表3 BP神經網絡模型中的變量

表4 韓家峁、白家川、趙石窯各序列擬合度R2對比

4 結論

在以往對水文時間序列自相關性的研究中,多聚焦于序列的一階自相關性,認為一階自相關成分剔除后,高階自相關即隨之消除或不再顯著。然而目前已有高階自相關的相關研究成果表明,其影響客觀存在且不容忽視。由于國內外類似研究較少,水文序列中的高階自相關性對水文分析是否有影響、有何種影響尚不明確,故本文對上述問題做了一些有益的嘗試與分析。具體為使用分段去趨勢方法對不同流域實測年徑流序列進行處理,并對去趨勢序列的高階自相關性進行識別;以BP神經網絡為基礎,降水-徑流模型為載體,探討了加入低階或(與)高階降雨序列對徑流預測效果的影響,得到了以下結論。

a) 高階自相關性不同于一階自相關性,其并不是一個普遍問題。不同流域徑流序列中高階自相關性的存在與否以及結構特征均有所差異。如本次研究中,漢江流域的石泉和安康站的實測年徑流序列無高階自相關性;渭河流域7個水文站的實測年徑流序列亦無高階自相關性;而無定河流域的韓家峁、橫山、趙石窯、白家川4個水文站的實測年徑流序列被檢測到高階自相關成分,其余7站則未發現高階自相關性。

b) 徑流高階自相關性是否存在,受地貌類型、土壤特性、降水規律、地下水循環條件等多種因素的共同影響。被檢測到高階自相關性的無定河流域的韓家峁、橫山、趙石窯水文站具有類似的影響特征,如:①均位于風沙區;②地下水補給徑流所占比重較大,一般達80%~90%以上;③陜北地區的包氣帶厚度一般較大,地下水由入滲到補給河流這一過程所需的時間較長,因此在實測徑流序列中呈現出了高階自相關性。漢江與渭河流域河川徑流主要來源于地表水,地下水基流比重較小,因而相關特征在水文站的實測年徑流序列表現不強烈,即不能檢測到顯著高階自相關性。

c) 水文序列中存在高階自相關性可能會影響統計檢測結果的準確性,進而得到不科學的結論。但高價自相關性并非只有劣勢,在徑流預測中如考慮水文序列高階自相關性影響,實施預測模型輸入變量優化,其可有效提高徑流預測預報的精度和效果。除此之外,在僅有短期實測徑流資料的情況下,高階自相關性的有效利用,也可為展延年徑流資料提供一種新思路。

猜你喜歡
趨勢
趨勢
第一財經(2025年5期)2025-05-16 00:00:00
退休的未來趨勢
英語世界(2023年12期)2023-12-28 03:36:16
趨勢
第一財經(2021年6期)2021-06-10 13:19:08
趨勢
初秋唇妝趨勢
Coco薇(2017年9期)2017-09-07 21:23:49
未來直銷的七大趨勢
趨勢
流行色(2016年10期)2016-12-05 02:27:24
SPINEXPO?2017春夏流行趨勢
關注醫改新趨勢
中國衛生(2015年2期)2015-11-12 13:14:02
“去編”大趨勢
中國衛生(2015年7期)2015-11-08 11:09:38
主站蜘蛛池模板: 日韩国产亚洲一区二区在线观看| 亚洲第一区欧美国产综合| 久久青草免费91线频观看不卡| 香蕉99国内自产自拍视频| 伊人国产无码高清视频| 乱人伦99久久| 久久综合亚洲鲁鲁九月天| 国产sm重味一区二区三区| 97影院午夜在线观看视频| 亚洲一级毛片免费看| 国产va免费精品观看| 亚洲中文字幕97久久精品少妇| a级毛片免费播放| 91精品国产91久无码网站| 中文字幕在线看视频一区二区三区| 精品在线免费播放| 午夜视频www| 欧美性久久久久| 精品久久久久无码| 国产乱子伦无码精品小说| 天天干天天色综合网| 成人无码一区二区三区视频在线观看| 一区二区三区国产精品视频| 重口调教一区二区视频| 成年人午夜免费视频| 欧美一区二区啪啪| 亚洲一级毛片在线观播放| 手机精品福利在线观看| 国产精品女熟高潮视频| 无码精品国产dvd在线观看9久| 久久国产高清视频| 亚洲不卡影院| 华人在线亚洲欧美精品| 亚洲精品国产日韩无码AV永久免费网 | 欧美成人h精品网站| 久久这里只精品国产99热8| 日本a级免费| 国产精品青青| www.日韩三级| 成人综合久久综合| 亚洲色成人www在线观看| 久久福利网| 99er精品视频| 日本欧美精品| 91精品啪在线观看国产| 色婷婷丁香| 自拍亚洲欧美精品| 亚洲人成人伊人成综合网无码| 久草视频福利在线观看| 国产呦精品一区二区三区网站| 国产成人精品亚洲77美色| 免费xxxxx在线观看网站| 亚洲人成网18禁| 国产精品片在线观看手机版 | 毛片a级毛片免费观看免下载| 欧美精品综合视频一区二区| 亚洲欧美国产五月天综合| 粗大猛烈进出高潮视频无码| 成人精品亚洲| 黄色成年视频| 久久青草免费91观看| 亚洲一级毛片在线观| 国产在线无码一区二区三区| 国产一级毛片在线| 国产精品视频观看裸模 | 亚洲色图欧美激情| 亚洲精品视频免费观看| 亚洲国内精品自在自线官| 97视频在线精品国自产拍| 欧美a级完整在线观看| 久草热视频在线| 欧美第九页| 99免费在线观看视频| 国产成熟女人性满足视频| 一本大道香蕉久中文在线播放| 波多野结衣二区| 亚洲一本大道在线| 国产精品3p视频| 欧美翘臀一区二区三区| 亚洲色中色| 尤物亚洲最大AV无码网站| 午夜日b视频|