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

黃河上游徑流豐枯空間分布特征及其影響因素

2021-12-06 09:50:22魏伊寧李勛貴
水資源保護 2021年6期
關鍵詞:模態影響

魏伊寧,李勛貴,李 芳

(1.蘭州大學資源環境學院,甘肅 蘭州 730000; 2.江西師大附屬外國語學校,江西 南昌 330103;3. 蘭州大學冰川與沙漠研究中心,甘肅 蘭州 730000; 4.廣西大學土木建筑工程學院,廣西 南寧 530004)

水分盈缺對陸地生態系統影響巨大,由水分盈缺引起的旱澇災害是目前最嚴重的災害之一[1]。進入21世紀后,中國年平均旱災與20世紀相比有不同程度的增加[2]。水文干旱是氣象干旱和農業干旱的延續與發展,因此被認為是最徹底的干旱[3-4],它反映了地表水與地下水的異常偏少現象,是受氣候變化和人類活動共同影響下的結果。水文干旱指標的構建通常基于徑流量,以時段徑流量小于某臨界值來表征[5]。標準徑流指數(standard runoff index,SRI)是指示水文干旱的重要指標,其本質反映了徑流豐枯的變化程度[6]。

黃河上游跨越青海、甘肅、陜西和內蒙古4個省、自治區,是南水北調西線工程得以實現和建設“新絲綢之路經濟帶”國家戰略的核心通道[7-8],具有特殊的地理環境和戰略意義。其中,蘭州以上黃河流域位于青藏高原東部,為生態環境脆弱區。水資源的高強度開發利用,改變了黃河上游徑流豐枯變化程度[9-10]。20世紀90年代,黃河上游連續極端枯水年的出現,加之流域水資源開發強度的持續增加,流域中下游斷流現象異常嚴重[11-13]。在全球變暖與人類活動的雙重影響下,徑流豐枯變化更加復雜,給區域日益嚴峻的水資源供需矛盾增加了新的危機[14-15]。但目前對于青藏高原背景下的區域徑流豐枯空間分布特征與其影響因素之間的聯系仍然不清楚,有很多問題亟待解決。因此,研究黃河上游徑流豐枯變化特征及其影響因素,對揭示區域水資源變化規律和保障區域水安全具有重要意義。

目前黃河流域徑流豐枯特征的研究重點側重于對時間、頻率特征的分析[16-18],而對徑流豐枯的空間分布特征研究較少。同時,已有研究多側重于采用某一水文站或多個水文站的平均情況來代表整個黃河流域,難以反映區域內部徑流豐枯的變化特征[19-20]。SRI能較為簡便有效地反映不同時間尺度下徑流的豐枯變化特征[21-23],目前在黃河上游應用較少。因此,本研究基于徑流SRI序列,結合經驗正交函數(empirical orthogonal functions,EOF)[24-26]將區域徑流豐枯的時間特征與空間結構分離,獲得有意義的空間分布模態來表征區域徑流豐枯的空間分布特征,然后運用奇異值分解(singular value decomposition,SVD)方法[27]揭示區域徑流SRI豐枯序列與存在時域相關性的青藏高原不同氣象因子場(降水、氣溫和積雪)之間的空間聯系。

1 研究區概況與數據來源

1.1 研究區概況

黃河是中國第二長河流,是中國西北和華北地區灌溉農業、生活、工業和生態服務的主要水源,但又一直面臨著水資源短缺和相關的環境問題[28-29]。研究區為蘭州以上的黃河上游地區,該區是黃河的主要產水區,也是中國十三大水電能源基地之一,流域面積為22.3萬km2,占流域總面積的29.6%,年徑流深度為147.1 mm,年平均流量為3.274×1010m3,占河流總流量的56.4%[30-31]。研究區集水面積大,跨度長,主要山脈海拔高程在2 800~4 800 m,山勢高聳,峽谷眾多,處于青藏高原和黃土高原的過渡區,海拔相差懸殊,地理環境復雜。研究區處于季風影響邊緣區,加之青藏高原的隆起對大氣環流有顯著影響,氣候要素的年、季變化大[32]。研究區濕度小、蒸發大,年平均降水量為426.2 mm,年蒸發量超過2 500.0 mm[33]。受梯級水庫及下游灌區影響,黃河流域的水文循環發生了很大變化,使得蘭州以上的黃河上游成為徑流豐枯變化最大、空間分布最不均勻的地區之一。

1.2 數據來源

本研究使用的數據為黃河上游蘭州以上地區7個典型水文站1969—2012年實測月徑流數據及青藏高原地區月降水量、月平均氣溫和月累積積雪深度數據。水文站位置如圖1(a)所示,分別為吉邁(JM)、瑪曲(MQ)、唐乃亥(TNH)、貴德(GD)、循化(XH)、小川(XC)和蘭州(LZ)。雪量站分布在10個流域,從西至東YDH、TLM、IN、YLZBJ、CDM、CJ、LCJ、NJ、HXZL和YH分別代表印度河、塔里木、內陸河、雅魯藏布江、柴達木、長江、瀾滄江、怒江、河西走廊和黃河流域(圖1(b))。實測徑流數據來源于中國黃河水利委員會(http://www.yrcc.gov.cn/);降水、氣溫和積雪深度數據分別來源于國家青藏高原數據中心(http://westdc.westgis.ac.cn/zh-hans/)提供的中國降水、氣溫格點數據集V2.0和青藏高原臺站雪深數據集V1.0。所用數據經交叉驗證和誤差分析,質量良好。

(a) 水文站分布

(b) 雪量站分布

對數據預處理時選取1969—2012年月徑流序列構建SRI序列,時間尺度采用季節(3個月)和年(12個月)尺度。對青藏高原月降水量、月平均氣溫和月累積積雪深度也先提取出不同時間尺度下(3個月、12個月)的時間序列,剪裁到研究區范圍內。將不同時間尺度下SRI、降水量、平均氣溫及積雪深度序列作為分析資料進行距平化處理后再標準化,以標準化后的距平值為對象進行EOF和SVD分析。

2 研究方法

2.1 SRI計算

SRI是由Shukla和Wood于2008年提出來的一種評價流域水文干旱和徑流豐枯程度的計算方法,步驟[18]如下:

步驟1:根據一定時間尺度的徑流分布概率密度函數fY(SRI指數的構建一般采用Gamma分布[21,34]),對其求累積概率F(Y):

(1)

式中Y為給定的時間尺度的徑流量。

步驟2:采用等概率轉換的方法,把累積分布概率F(Y)轉換為均值為0、方差為1的正態分布,即得到一定時間尺度范圍的SRI值,計算公式為

(2)

其中

I=[lnF(Y)-2]0.5

式中:VSRI為SRI值,VSRI越小表示徑流越小(枯水)或干旱,反之則表示徑流越大(豐水)或濕潤;參數設置參考趙茹欣等[35]的研究,取c0=2.515 517,c1=0.802 853,c2=0.010 328,d1=1.432 788,d2=0.189 269,d3=0.001 308。

2.2 EOF方法

對得到的SRI序列進行EOF分解,步驟[36]如下:

步驟1把得到的SRI值預處理后排列為矩陣X(有m個空間點和n個時間點數),計算X與轉置矩陣XT的交叉積,得到方陣C,對C進行EOF分解,得到特征根矩陣λ和分解矩陣W。

步驟2定義λj為第j模態的特征根,得到第j模態的EOF值Wj:

Wj=(W1j,W2j,…,Wmj)T

(3)

步驟3將Wj投影到原始資料矩陣X上,就得到所有空間特征向量對應的時間系數矩陣P:

(4)

2.3 SVD方法

Prohaska最先將SVD應用于氣候學研究,用于診斷美國月平均地面氣溫與北太平洋海平面氣壓之間的關系[37],其計算過程如下:

步驟1計算預處理后觀測矩陣X(有m個空間點和n個時間點數)與氣象因子矩陣Y(有q個空間點和n個時間點數)的交叉協方差陣Q,對其進行奇異值分解,得到U、V空間奇異向量。觀測矩陣X也稱左場,氣象因子矩陣Y也稱右場。

步驟2設第j模態空間奇異向量為Uj和Vj。把原觀測場X投影到Uj,把Y場投影到Vj,即可得到第j模態的左右場的時間序列,記為aj和bj:

(5)

步驟3計算X、Y場每個空間觀測點的時間序列(X、Y的某一行)分別與時間系數aj、bj的異性相關系數r:

(6)

式中:ajN、bjN分別表示第j模態第N個時間點的時間系數;σj為特征根矩陣Σ第j模態的特征值。

2.4 方差貢獻率與顯著性檢驗

a.用方差貢獻率表征第j模態解釋的協方差占總協方差的比例,記為QSCF:

(7)

式中:M為總模態數;Qj為第j模態特征值;Qk為第k模態特征值。

b.為避免出現“虛假模態”,本文采用North檢驗對EOF結果進行顯著性檢驗:

(8)

當λj-λj+1≥ej時,認為EOF分解出的模態為獨立的模態。

c.為避免出現“虛假相關”,本文采用t檢驗對SVD相關系數結果進行顯著性檢驗:

(9)

式中:α為置信度,本研究取為95%;tα為置信度為α的t分布的上α分位數。如果相關系數r的絕對值大于置信度為95%的Rc,則認為通過顯著性檢驗。

3 結果與分析

3.1 SRI的空間分布

3.1.1年尺度

式(3)中將第1個特征值對應的特征向量在空間上的分布稱作第1空間分布模態,即W1;第2個特征值對應的特征向量在空間上的分布稱為第2空間分布模態,即W2;以此類推。特征向量(即EOF值)的正負表示正反位相,特征向量值越大認為這個空間點距離平均值越大,即越敏感。計算結果如圖2所示,可見黃河上游地區徑流豐枯變化有兩個典型的分布模態,均通過了North檢驗,檢驗結果見表1。

(a) 第1空間分布模態

(b) 第2空間分布模態

表1 SRI的方差貢獻率和North檢驗

由圖2和表1可見,年尺度徑流SRI的第1模態的方差貢獻率達到85.2%,為黃河上游徑流SRI的最主要空間分布模態。由圖2(a)可見,整個研究區的EOF都為正值,說明整個區域徑流SRI在一定程度上受到某些共同因子的影響而表現出一致的主要變化趨勢,即某些年份區域徑流量整體趨于偏豐或偏枯狀態。由圖3可見,第1空間分布模態的時間系數在1989年以前大多數年份為正值,表示全區水量偏豐,且在1975年和1983年表現最突出;1989—2002年間大多數年份為負值,表示全區水量偏枯,且在1997年和2002年是表現最突出,說明黃河上游水文干旱最為嚴峻的時間段為20世紀90年代,與現實情況相符。由圖2(b)可見,第2空間分布模態的方差貢獻為10.7%,為黃河上游徑流SRI的次要空間分布模態。該模態下徑流SRI序列呈東北—西南的反位相分布特征,即某些年份流域東北部徑流量偏豐,而西南部徑流量則偏枯,反之亦然。由圖3可見,第2空間分布模態時間系數在1986年以前在0值附近波動,說明1986年以前這一空間分布模態并不典型,主要以第1空間分布模態(即全區變化具有一致性)為主;1986年后時間系數呈正負交替振蕩,振蕩幅度加大,說明1986年后第2空間分布模態開始變得顯著,且流域東北、西南部徑流豐枯轉變較為頻繁。

圖3 年尺度SRI空間分布模態的時間系數

(a) 春季第1空間分布模態

(b) 夏季第1空間分布模態

(c) 秋季第1空間分布模態

(d) 冬季第1空間分布模態

(e) 春季第2空間分布模態

(f) 夏季第2空間分布模態

(g) 秋季第2空間分布模態

(h) 冬季第2空間分布模態

(i) 冬季第3空間分布模態

(b) 第2空間分布模態

3.1.2季節尺度

季節尺度SRI序列除冬季有3個顯著的分布模態(圖4(i))外其他季節均具有2個顯著的分布模態。各季節的空間分布特征如圖4所示,圖中顏色更深的區域表示更敏感。

由圖4可見,季節尺度SRI的空間分布整體上與年尺度的空間分布較為相似,4個季節的第1空間分布模態均呈現出全區同位相變化特征,第2空間分布模態均呈現出東北—西南反位相變化特征,但各季節敏感區有所不同。由第1空間分布模態的時間系數來看(圖5(a)),冬季和春季均于2003年負值最為顯著,說明2003年冬、春季黃河上游徑流水文干旱情勢最為嚴峻;夏季和秋季具有明顯的時間變化的階段性特征,說明黃河上游夏、秋季可能出現多年偏豐或偏枯的情形,且時間系數呈現出上升(偏豐)—下降(偏枯)—上升(偏豐)的趨勢,下降最強烈的時期為20世紀90年代,說明20世紀90年代區域徑流持續偏枯,與實際情況相符。由第2空間分布模態時間系數來看(圖5(b)),4個季節都具有較明顯的階段特征,冬季和春季變化較為類似,呈波動上升趨勢,且近年來以正值為主,說明近年來SRI的冬、春季空間分布的次要模態呈現出西南區域水量偏枯而東北區域水量偏豐的特征;夏季和秋季的時間系數變化呈波動下降趨勢,且近年來以負值為主,說明近年來SRI的夏、秋季空間分布次要模態呈現出西南區域水量偏豐,東北區域水量偏枯的特征。

3.2 SRI與氣象因子耦合空間分布特征

3.2.1年尺度

由式(6)可計算SRI序列與青藏高原氣象因子的異性相關系數。SRI與降水、氣溫、積雪的相關系數分別為0.85、0.40、0.47,并都通過置信度95%的顯著性檢驗。整體來看,黃河上游地區SRI與降水場有最好的相關關系,積雪次之,氣溫最弱。除少數幾個季節未通過顯著性檢驗外,其他均通過置信度為95%的t檢驗,說明青藏高原降水、積雪和氣溫均是影響黃河上游徑流量豐枯變化的關鍵氣象因子。

不同氣象因子影響流域的關鍵區不同,圖6和圖7為年尺度SRI與不同氣象因子相關系數空間分布與時間系數。可見,青藏高原東北部流域的降水(圖6(d))、北部氣溫(圖6(e))以及東部積雪(圖6(f))對整個黃河上游流域徑流豐枯變化有顯著影響。黃河上游SRI與各氣象因子的空間聯系均表現出一致的變化規律(圖6(a)(b)(c)),與SRI的第1模態空間分布類似,說明研究區SRI的一致偏豐或偏枯的空間分布模態是受氣象因子的影響形成的。

3.2.2季節尺度

從表2可以看出,不同季節SRI與氣象因子耦合場的相關性不同,相關性最大的為夏季降水與夏季SRI,相關系數達0.81,夏季降水與后期秋季SRI的相關系數也較大,達0.71,說明夏季降水除了對同期SRI有影響較大外,還具有一定的時間滯后性。降水的季節尺度相關系數均沒有年尺度的相關系數大,說明降水對SRI的影響持續時間較長。從氣溫方面來看,夏季氣溫與春季SRI相關系數最大,達0.70,且大于年尺度相關系數。從積雪方面來看,春季積雪與春季SRI的相關系數最大,達0.64,且大于年尺度相關系數。氣溫、積雪對SRI季節性影響較年尺度明顯,表明氣溫、積雪對SRI的影響持續時間較短。

表2 不同氣象因子與SRI的相關系數

選取季節尺度的空間分布相關系數較大的季節進行分析(圖8),可以發現降水、氣溫對徑流的主要影響區都有擴大的趨勢(圖8(d)(e)),積雪的主要影響區(圖8(f))由東側轉向青藏高原中部,即黃河、長江流域上游靠上位置,說明積雪作為青藏高原下墊面的一種強迫源,季節性消漲較大,主要通過影響黃河源頭季節性(春季)積雪融水來影響上游徑流的豐枯的分布。

4 討 論

4.1 SRI與其他干旱指標的對比

目前,利用干旱指數分析黃河流域徑流豐枯變化的研究中,大多基于氣象干旱指數,從降水、蒸發角度研究流域徑流豐枯狀況,如周帥等[38]基于標準化降水指數表明夏季黃河源區旱情嚴重,干旱歷時長、烈度大;牛亞婷等[39]采用標準化降水指數揭示了秋季黃河流域西部多發生特旱,與本研究夏季、秋季區域更容易出現多年徑流量偏豐或偏枯的情形類似;WANG等[40]和HUANG等[41]均指出黃河流域1990年代干旱災害尤為嚴重,季節干旱化趨勢明顯,與本研究的20世紀90年代區域一致的徑流量偏枯結論一致;佘敦先等[42]在識別出黃河流域極端干旱事件的基礎上,采用Copula函數構建了干旱歷時和干旱強度兩變量統計模型;劉勤等[43]采用相對濕潤度指數指出黃河上游旱情最為嚴重,由西北至東南方向旱情逐漸減弱。但這些研究所采用的干旱指數均基于單站點,忽視了不同區域的空間變化特征。由于徑流-降水具有一定的時間滯后性,形成枯水的時間也會具有一定的時間延遲,從氣象角度則難以反映出來。因此,從徑流的角度建立的指標能更直接地反映徑流豐枯變化特征。現階段,從徑流角度構建的指標有距平百分率、Z指數、濕度指標等。距平百分率方法簡單,但會出現相同的指標值對應不同的豐枯狀況的現象;Z指數對豐枯情況響應快,但小區域的適用性不強;濕度指標相對上述兩種指標而言過于敏感,因而可能存在過分夸大了豐枯程度的情況[44]。SRI指標計算簡單且能較為準確地反映流域徑流豐枯狀況;EOF和SVD方法則是研究地理要素空間分布特征的常用方法。因此,選用SRI與EOF、SVD耦合方法能較為準確地描述黃河上游徑流豐枯的空間分布特征及與氣象因子的空間聯系。

(a) SRI與降水左場

(b) SRI與氣溫左場

(c) SRI與積雪左場

(d) SRI與降水右場

(e) SRI與氣溫右場

(f) SRI與積雪右場

(a) 左場時間系數

(b) 右場時間系數

(a) 夏季SRI與夏季降水左場

(b) 春季SRI與夏季氣溫左場

(c) 春季SRI與春季積雪左場

(d) 夏季SRI與夏季降水右場

(e) 春季SRI與夏季氣溫右場

(f) 春季SRI與春季積雪右場

4.2 影響SRI空間分布的原因探討

4.2.1氣象因素

氣象因素對黃河上游徑流豐枯的空間分布有重要影響,降水是影響徑流豐枯分布最關鍵的氣象因子。除降水外,氣溫和積雪與SRI的季節最大相關系數均大于年尺度的相關系數,說明氣溫和積雪對徑流豐枯分布的影響時間尺度較短,僅在某個季節具有顯著的影響;而降水對徑流豐枯分布的影響既有短期的時間尺度上的顯著影響,又由于其具有一定的時滯效應和累積效應使其在較長時間尺度上也有較大作用。結合SRI分布特征及其與氣象因子的關系,可探究影響黃河上游徑流豐枯分布特征的氣象驅動機制。在全球氣候變暖背景下[45],青藏高原年平均溫度從20世紀80年代中期開始呈持續上升趨勢[46],與本研究圖7(b)、圖9(b)時間序列趨勢一致。青藏高原北部由于受到緯度及反照率的影響對溫度變化響應更加敏感(圖6(e)),與徐國昌等[47]研究結果一致,該區域氣溫的變化不僅影響了黃河流域源頭春季融雪水量的變化,同時積雪異常會加劇大氣的熱力狀況的變化程度,從而影響該地區大氣環流狀況[48-50]。在氣溫較高年份,積雪對大氣的冷卻作用減弱導致了高原大氣熱源進一步加強,這不僅引導了印度季風向高原內部的延伸,也使得東亞夏季風越過高大山脈及地形阻擋,給青藏高原東北側的黃河上游帶來更多的降水,出現區域一致性水量偏豐的情況[51]。而一般年份,來自太平洋的濕潤水汽難以到達青藏高原內陸和北部,且青藏高原東北部位于氣流上升運動的下沉補償區[52-54],下沉氣流較干燥,此時區域降水較少,故黃河上游出現一致性水量偏枯的情況。

(a) 左場時間系數

(b) 右場時間系數

4.2.2人類活動

蘭州以上黃河流域的人類活動影響主要在于梯級水庫的運用,尤其是多年調節水庫——龍羊峽水庫和年調節水庫——劉家峽水庫的影響。唐乃亥站上游由于人類活動影響較小,其實測徑流可近似地認為是天然徑流。唐乃亥站和貴德站之間沒有大的支流匯入,且人類引用水也較少,故理論上如果沒有龍羊峽水庫的調節,唐乃亥站和貴德站的空間分布模態應該相似,而圖2(b)結果顯示,黃河流域上游具有兩個變化中心,即唐乃亥—瑪曲—吉邁和貴德—循化—小川—蘭州,唐乃亥站和貴德站的徑流變化特征相差較大,兩者的差異可以認為是龍羊峽水庫調節引起的。且這一空間分布模態在1986年后表現更加顯著,正好與1986年龍羊峽水庫蓄水時間點一致,使得1986年后出現東北—西南兩極分化的反位相空間分布模態,流域東北、西南區徑流豐枯轉變也變得頻繁。這表明人類活動對黃河上游徑流豐枯的空間分布存在影響。

5 結 論

a.黃河上游徑流豐枯的空間分布有兩種典型的分布模態,一種是全區一致性分布模態,一種是東北—西南兩級分化的反位相分布模態。以一致性分布模態為主,在20世紀90年代最典型,且夏季、秋季更容易出現全流域的多年徑流量偏豐或偏枯的情形。

b.黃河上游徑流豐枯的空間分布是受多因素影響的復雜過程,氣象因子和人類活動的綜合作用是其變化的根本成因。一致性分布模態主要受氣象因子的影響,兩級化分布模態受人類活動的影響較大。

c.影響黃河上游徑流豐枯變化的氣象因子中,降水是主要因子,降水的變化對徑流豐枯變化具有較大的持續性影響,而氣溫、積雪的季節性影響較明顯。影響黃河上游徑流豐枯變化的人類活動因子中,水庫的調節是主要因子。

d.不同氣象因子影響黃河上游徑流豐枯變化的關鍵區不同,青藏高原東北部的降水、東部的積雪及北部的氣溫對黃河上游徑流豐枯變化具有顯著影響,且春、夏季降水和氣溫的關鍵區會擴大,積雪的關鍵區會轉移。

猜你喜歡
模態影響
是什么影響了滑動摩擦力的大小
哪些顧慮影響擔當?
當代陜西(2021年2期)2021-03-29 07:41:24
沒錯,痛經有時也會影響懷孕
媽媽寶寶(2017年3期)2017-02-21 01:22:28
車輛CAE分析中自由模態和約束模態的應用與對比
擴鏈劑聯用對PETG擴鏈反應與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
基于Simulink的跟蹤干擾對跳頻通信的影響
國內多模態教學研究回顧與展望
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 91色在线观看| 中国一级特黄视频| 久久伊伊香蕉综合精品| 毛片a级毛片免费观看免下载| 中文字幕第1页在线播| 在线观看视频一区二区| www.91中文字幕| 久久精品国产精品青草app| 国产人前露出系列视频| 国产超碰在线观看| 欧美在线精品怡红院| 91九色国产porny| 黄色网站在线观看无码| 国产国产人免费视频成18| 国产在线观看成人91| 91蝌蚪视频在线观看| 亚洲αv毛片| 人妻精品全国免费视频| 国产熟女一级毛片| 国产一二视频| 国产成人乱无码视频| 高清不卡一区二区三区香蕉| 欧美色视频日本| 亚洲欧美成人综合| 91小视频在线观看| 亚洲成人黄色网址| 一级在线毛片| 亚洲天堂久久| 国产成人综合在线视频| 国产91高清视频| 女人18毛片一级毛片在线 | 在线人成精品免费视频| 午夜国产精品视频黄| 沈阳少妇高潮在线| 极品国产在线| 国产一区二区人大臿蕉香蕉| 黄片在线永久| 亚洲天堂日本| 亚洲第一成人在线| 无码中字出轨中文人妻中文中| 欧美国产中文| 一级毛片在线免费看| 国产精品久久久久久久久kt| 亚洲成av人无码综合在线观看| 国内99精品激情视频精品| 国产精品亚洲欧美日韩久久| 欧美日本在线观看| 国产在线自在拍91精品黑人| 无码有码中文字幕| 国产第一色| 国产欧美日韩在线一区| 亚洲视频免费播放| 亚洲综合九九| 国产午夜不卡| 国产特级毛片| 六月婷婷综合| 67194亚洲无码| 国产精品免费入口视频| 亚洲精品无码抽插日韩| 一区二区三区精品视频在线观看| 一级不卡毛片| 成人午夜网址| 亚洲欧洲自拍拍偷午夜色| 久久国产成人精品国产成人亚洲| 2021精品国产自在现线看| 精品伊人久久久久7777人| 色综合久久无码网| h视频在线播放| 国产伦片中文免费观看| 热99精品视频| 亚洲日韩精品欧美中文字幕 | 白浆免费视频国产精品视频| 青草精品视频| 996免费视频国产在线播放| 日韩视频福利| 国产成人a在线观看视频| 天堂成人av| 一区二区三区成人| 国产午夜不卡| 欧美国产菊爆免费观看| 91在线播放免费不卡无毒| a级毛片网|