• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    黃河寧夏段干支流非一致性洪峰遭遇風(fēng)險分析

    2016-12-08 08:24:32李子遠苑希民
    水利水電科技進展 2016年6期
    關(guān)鍵詞:兩河清水河洪峰

    李子遠,馮 平,苑希民

    (天津大學(xué)水利工程仿真與安全國家重點實驗室,天津 300072)

    ?

    黃河寧夏段干支流非一致性洪峰遭遇風(fēng)險分析

    李子遠,馮 平,苑希民

    (天津大學(xué)水利工程仿真與安全國家重點實驗室,天津 300072)

    對黃河和清水河非一致性洪峰進行遭遇風(fēng)險分析。采用綜合診斷法對洪峰序列進行變異性診斷,引入混合von Mises 分布和混合P-Ⅲ分布對洪峰發(fā)生時間和洪峰量級進行擬合;選擇非對稱部分嵌套Archimedean Copula函數(shù),構(gòu)建兩河洪峰聯(lián)合分布并分析兩河洪峰發(fā)生時間和量級遭遇情況,繪制各種頻率組合遭遇概率分布圖。結(jié)果表明:兩河洪峰序列均發(fā)生變異,兩河洪峰發(fā)生時間遭遇概率在7月31日和8月31日達到最大,分別為3.36×10-4和1.15×10-4;不同頻率洪峰遭遇概率在7月26日至8月2日和8月31日至9月5日兩個時段達到最大。

    洪峰遭遇;非一致性;混合分布;Copula函數(shù);風(fēng)險分析;黃河;清水河

    特大洪水往往是由多來源洪水共同作用形成的,特別是流域性大洪水,各干支流洪水之間的相互影響尤為突出。目前水文工作中常用方法是將各來源洪水視為不相關(guān),對洪水實測資料進行統(tǒng)計分析,該做法忽略了同流域內(nèi)不同來源洪水之間的內(nèi)在聯(lián)系,所得結(jié)果顯然不能準確描述實際洪水特征,且只能定性分析干支流遭遇情況,無法分析設(shè)計洪水情況。因此,綜合考慮干支流洪水內(nèi)在聯(lián)系,定量分析干支流洪水遭遇概率,對河道防洪和水庫調(diào)度有著重要意義。

    Copula函數(shù)能很好地描述多變量之間的相依性,作為連接函數(shù),具有形式多樣化且不受分布函數(shù)限制、在變量單調(diào)變換下形式不發(fā)生變化等特點,被越來越廣泛地用來對存在相關(guān)性的洪水序列進行聯(lián)合分析[1]。Yan等[2]通過構(gòu)建非對稱三變量Clayton Copula函數(shù)對南水北調(diào)中線供水區(qū)和受水區(qū)降水豐枯遭遇概率進行分析,計算了極端降水條件下輸水概率;Genest等[3]構(gòu)建Meta-elliptical Copula函數(shù)對Romaine河春汛的洪峰、洪量和持續(xù)時間進行了聯(lián)合頻率分析;Chen等[4]構(gòu)建X-Gumbel Copula函數(shù)分別對長江上游和Colorado河上游洪峰遭遇概率進行了分析。

    氣候環(huán)境變化和人類活動的影響導(dǎo)致部分流域水文序列的一致性遭到破壞,直接用實測數(shù)據(jù)進行邊緣分布擬合會降低計算結(jié)果可靠性,因此必須進行水文序列變異性檢驗,并對發(fā)生變異的序列進行一致性處理。Zhang等[5]采用Mann-Whitney U檢驗法和累積離差系數(shù)分別對珠江流域的西江、北江和東江年最大洪峰序列進行了變異診斷,并分析了變異點前后洪峰序列變化趨勢。目前常從兩個角度進行水文序列的非一致性分析,一是基于還原/還現(xiàn)將水文序列修正到過去狀態(tài)或修正到現(xiàn)在狀態(tài),以保證序列的一致性;二是直接對水文序列進行頻率分析[6]。

    洪水遭遇可定義為各來源洪水的最顯著特征量間的聯(lián)合遭遇(如洪峰、洪量和持續(xù)時間等)[7],相對于分析實測洪水序列的實際遭遇情況,以顯著特征量表征洪水序列在抓住所研究洪水特性的同時,能更靈活有效地使用聯(lián)合分析方法。

    本文以黃河和清水河為研究對象,選擇洪峰流量為洪水顯著特征量進行干支流洪水遭遇分析。對干支流洪水序列進行變異性診斷,采用混合分布擬合非一致性洪水序列,擬合年最大洪峰發(fā)生時間和量級分布;引入Archimedean Copula函數(shù),構(gòu)建黃河和清水河年最大洪峰發(fā)生時間和量級的聯(lián)合分布,進而分析兩河非一致性洪峰在發(fā)生時間和量級上的遭遇概率。

    1 研究區(qū)概況

    黃河寧夏段自中衛(wèi)南長灘入境,流經(jīng)衛(wèi)寧灌區(qū)、青銅峽水庫、青銅峽灌區(qū),至石嘴山麻黃溝出境,區(qū)內(nèi)流程397 km,分為峽谷段、庫區(qū)段和平原段,區(qū)域水系分布如圖1所示。入境洪水來自黃河上游流域降水。20世紀以來,發(fā)生洪峰流量超過6 000 m3/s的大洪水4次,分別為1904年、1946年、1964年和1981年,其中1904年和1946年對寧夏沿黃兩岸造成嚴重洪災(zāi),而1964年和1981年由于預(yù)報及時和采取了有效的防洪措施,很大程度上減少了洪災(zāi)損失。區(qū)間主要有清水河、紅柳溝和苦水河3條支流匯入,其中對干流流量影響較大的是清水河,其流域面積達14 481 km2。清水河水系降水量為200~600 mm,由南向北遞減;降水的年際變化大,最大相差3倍,年內(nèi)分配也不均勻,多集中在7—9月。清水河多發(fā)暴雨洪水,且以局部暴雨為主,洪水過程一般持續(xù)3天,峰高量小,發(fā)生較大洪水的年份有1926年、1933年、1964年、1992年,其中1926年和1933年對下游造成較嚴重的洪水災(zāi)害,而1964年和1992年由于多座水庫攔洪,使得下游并未發(fā)生洪災(zāi)。在實測資料中,黃河和清水河并未發(fā)生大洪水遭遇(1964年二者均發(fā)生較大洪水,洪峰發(fā)生時間卻相差22 d),但不能因此忽略二者遭遇的可能性。目前青銅峽水庫淤積日趨嚴重,防洪能力減弱,使下游兩岸灌區(qū)的防洪形勢更加嚴峻,而黃河和清水河洪峰是否遭遇對青銅峽水庫入庫洪峰流量影響很大,因此對兩河洪峰遭遇情況進行分析顯得尤為重要。

    圖1 黃河寧夏段水系示意圖

    2 研究方法

    2.1 洪水序列非一致性分析

    目前國內(nèi)外對水文變異研究較多,Prohaska等[7-9]基于水文序列的確定性成分和隨機性成分,以水文序列分布形式或分布參數(shù)發(fā)生顯著變異為突變標準,提出了水文變異診斷的綜合診斷系統(tǒng)。

    a. 初步診斷。繪制洪水序列的過程線和5點滑動平均線,觀察序列的走勢,進而用R/S法和極差分析法計算洪水序列的Hurst系數(shù)h,并根據(jù)表1(表中hα=0.661 和hβ=0.703分別為α=0.05和β=0.01時h的臨界值)判斷洪水序列的強弱程度。

    表1 Hurst系數(shù)變異程度分級表[9]

    b. 詳細診斷。采用下列較為成熟的8種檢驗方法從不同角度對初步診斷判定存在變異的序列進行詳細診斷:MWP鑒定法[10]、滑動T檢驗法[11]、滑動F檢驗法、B-F檢驗法[12]、滑動秩和檢驗法[13]、L-H檢驗法、Cramer檢驗法[14]、有序聚類法[15]。由于不同檢驗方法假設(shè)前提和判斷標準不一致,得到的檢驗結(jié)果往往存在一定差異。

    c. 綜合診斷。詳細診斷是從數(shù)理統(tǒng)計角度找出可能的變異點,其結(jié)果不一定符合流域?qū)嶋H情況,還需結(jié)合物理成因分析,即下墊面變化、水利工程建設(shè)和水土保持工程等對水文序列的影響,對詳細診斷結(jié)果進行合理性甄別,最終確定變異點。

    2.2 洪峰分布形式確定

    由于洪峰遭遇分為洪峰發(fā)生時間遭遇和洪峰量級遭遇,將洪峰序列分為洪峰發(fā)生時間序列和洪峰流量序列,分別構(gòu)建洪峰發(fā)生時間分布和洪峰流量分布,基于變異診斷結(jié)果,對于存在變異的洪峰序列,采用多參數(shù)混合分布擬合。

    2.2.1 洪峰發(fā)生時間分布

    von Mises分布能很好地描述具有周期性或季節(jié)性的有限變量[16],選擇von Mises分布對洪峰發(fā)生時間進行擬合,以汛期起始和結(jié)束時間為上下邊界,將洪峰發(fā)生時間序列轉(zhuǎn)換為[0,2π]圓周上的方向序列,轉(zhuǎn)換公式為

    (1)

    式中:L為汛期長度;Di為洪峰發(fā)生時間。

    von Mises分布概率密度函數(shù)為

    (2)

    式中:μ為位置參數(shù),0≤μ≤2π;k為尺度參數(shù),k≥0;I0(k)為0階變形Bessel函數(shù)。

    von Mises分布為圓周上的單峰分布,而洪峰時間序列通常呈多峰形態(tài),單一分布函數(shù)往往不能準確描述序列分布形式,因此需要選擇多個von Mises分布進行混合分布擬合,形式如下:

    ,

    0≤μi≤2π,0≤ki)

    (3)

    式中:N為序列長度;αi為權(quán)重系數(shù),滿足α1+α2+…+αN=1。

    采用改進似然比檢驗(MLRT)[18]對方向數(shù)據(jù)進行檢驗,確定其分布形態(tài),構(gòu)造統(tǒng)計量:

    式中:n為數(shù)據(jù)個數(shù);R和θ分別為方向數(shù)據(jù)對應(yīng)的半徑和角度。當n較大時,近似服從χ2分布。

    由于von Mises分布的非正規(guī)性,通常情況矩估計不是漸進有效估計,極大似然估計得到結(jié)果并不是被估參數(shù)的相合估計,陳家驊等[17]對極大似然估計進行了改進。通過對比發(fā)現(xiàn)帶懲罰的極大似然估計的結(jié)果最優(yōu),因此本文基于改進的期望最大(EM)算法采用帶懲罰的極大似然估計進行混合von Mises分布參數(shù)估計。

    2.2.2 洪峰流量分布

    根據(jù)SL 44—2006《水利水電工程設(shè)計洪水計算規(guī)范》,選擇P-Ⅲ型分布擬合洪峰流量序列,其密度函數(shù)形式如下:

    (5)

    式中ε、β和b0分別為形狀、尺度和位置參數(shù)。

    對于存在變異的洪峰序列,單一分布形式已不能準確擬合實際序列。成靜清等[19]通過構(gòu)建兩個對數(shù)正態(tài)分布的混合分布對陜北實測平均流量序列進行擬合,得到比單一分布更好的成果。因此本文構(gòu)建多分布形式混合分布對兩河洪峰序列進行擬合,混合分布的密度函數(shù)形式為

    (6)

    式中:φi為權(quán)重系數(shù);εi、βi和b0i為各分項密度函數(shù)參數(shù)。對混合分布參數(shù),模擬退火算法是一種有效算法,然而傳統(tǒng)模擬退火算法存在一些局限,如計算步長難確定,收斂速度慢,易漏掉最優(yōu)解,為此眾多學(xué)者進行了改進,朱顥東等[25]對退火過程和抽樣過程進行了兩階段改進,并模擬得到比傳統(tǒng)算法更優(yōu)的解。因此,本文先用矩法進行參數(shù)初估計,再用模擬退火算法進行參數(shù)最終估計。

    2.3 聯(lián)合分布構(gòu)建

    Archimedean Copula函數(shù)族因其表達式簡單和易于構(gòu)造而被廣泛應(yīng)用,其中具有代表性的有Frank、Clayton和Gumbel函數(shù)。兩河洪峰序列具有上尾相關(guān)性,且洪峰遭遇屬于極值事件,Gumbel Copula函數(shù)能很好地滿足以上特征,因此本文選擇Gumbel Copula函數(shù)構(gòu)建聯(lián)合分布。兩變量聯(lián)合分布的表達式為

    C(u1,u2)=exp{-[(-lnu1)θ+(-lnu2)θ]1/θ}

    (7)

    式中:u1和u2為邊緣分布函數(shù);θ為連接參數(shù),θ>1。這里采用非參數(shù)法估計連接參數(shù)。

    進行洪水遭遇分析時,通常將不同河流的洪峰流量視為不相關(guān)[20]或?qū)⑼恿鞯暮榉灏l(fā)生時間和洪峰流量量級視為不相關(guān)[4],如此考慮固然能簡化計算,但忽略了洪峰序列之間的部分相關(guān)性,使得計算結(jié)果不能準確反映洪峰遭遇特征。本文在綜合考慮各河流洪峰序列自身內(nèi)在的和彼此間相關(guān)性的基礎(chǔ)上,構(gòu)建洪峰序列四變量聯(lián)合分布。由于每條河流洪峰發(fā)生時間序列和流量量級序列相關(guān)性較強,而河流之間的相關(guān)性較弱,相比于對稱Copula函數(shù),非對稱Copula函數(shù)能更好地描述多變量之間的相關(guān)性,且對數(shù)據(jù)的擬合更靈活[21]。選擇非對稱構(gòu)造中的半嵌套構(gòu)造,其表達式為

    C1(u1,u2)=exp{-[(-lnu1)θ1+(-lnu2)θ1]1/θ1}

    (8)

    C2(u3,u4)=exp{-[(-lnu3)θ2+(-lnu4)θ2]1/θ2}

    (9)

    C(u1,u2,u3,u4;θ1,θ2,θ3)=exp{-[(-lnC1)θ3+

    (-lnC2)θ3]1/θ3}

    (10)

    式中:u1、u2、u3、u4為邊緣分布函數(shù);θ1、θ2、θ3為連接參數(shù),且滿足1<θ3<θ1,1<θ3<θ2。采用偽極大似然估計法[24]進行參數(shù)估計,用樣本的經(jīng)驗分布函數(shù)代替邊緣分布函數(shù),則不需要進行邊緣分布函數(shù)參數(shù)估計,從而有效避免邊緣函數(shù)分布形式對相關(guān)參數(shù)估計結(jié)果的影響。

    采用Gringorten公式計算經(jīng)驗聯(lián)合頻率:

    (11)

    式中:m為聯(lián)合觀測樣本中滿足(xj

    選擇非參數(shù)K-S檢驗法,對聯(lián)合分布函數(shù)進行擬合檢驗,構(gòu)造統(tǒng)計量:

    (12)

    式中F0(x,y)為理論聯(lián)合分布概率。

    2.4 洪水遭遇風(fēng)險分析

    若兩河洪峰不在同一天發(fā)生,洪峰流量無遭遇可言,因而,需先對洪峰發(fā)生時間遭遇進行分析,并以此為中間量進行洪峰流量遭遇分析。洪峰發(fā)生時間遭遇風(fēng)險可定義為第i天同時發(fā)生洪水的概率:

    Pi=P(ti-1

    (13)

    式中T1和T2分別為兩河洪水發(fā)生時間。

    第i天兩河同時發(fā)生指定量級洪峰(q1,q2)的概率為

    P(i,q1,q2)=P(ti-1

    Q1>q1,Q2>q2)

    (14)

    式中Q1和Q2分別為兩河發(fā)生的洪峰流量。

    3 應(yīng)用分析

    3.1 干支流洪峰序列非一致性分析

    根據(jù)綜合診斷法的思路對黃河和清水河洪峰序列進行變異性診斷,找出序列中最大可能變異年份。通過水文序列的合理性分析,已將清水河泉眼山站洪峰序列1955年、1956年數(shù)據(jù)作為歷史特大值考慮,不參與變異性診斷。

    a. 初步診斷。黃河和清水河年最大洪峰過程線及5點滑動平均線如圖2所示。從圖2可知,黃河年最大洪峰的5點滑動平均線在1986年后呈下降趨勢;清水河年最大洪峰的5點滑動平均線在1993年前呈上升趨勢,1993年后呈下降趨勢。計算兩河Hurst系數(shù),黃河為0.858 7,屬于強變異;清水河為0.807 7,屬于中變異。初步診斷兩河洪峰時間序列存在變異。

    圖2 黃河和清水河洪峰序列及5點滑動平均

    b. 詳細診斷。采用選定的8種檢驗方法對黃河和清水河洪峰序列進行檢驗,最有可能變異點見表2。由表2可知,每種檢驗法均檢驗出1986年為黃河序列可能變異點且在各可能變異點中占權(quán)重最大,1991年為清水河序列可能變異點且在各可能變異點中占權(quán)重最大,可以斷定黃河在1986年發(fā)生變異,清水河在1991年發(fā)生變異。計算得黃河1986年前后和清水河1991年前后子序列的Hurst系數(shù)分別為0.527 4、0.541 3和0.616 2、0.587 3,均滿足一致性要求。

    表2 黃河和清水河變異點診斷結(jié)果

    c. 綜合診斷。黃河上游1986年龍羊峽水庫和劉家峽水庫實行聯(lián)合調(diào)度,對上游天然徑流進行攔蓄重分配,下游徑流過程發(fā)生很大變化,洪峰流量削減尤為突出;20世紀80年代末至90年代初,對清水河流域水利工程進行清淤整治,使得水庫群對洪峰的攔蓄調(diào)節(jié)能力顯著增強,特別是對流域內(nèi)經(jīng)常發(fā)生的暴雨洪水有顯著控制作用;20世紀80年代以來,隨著大氣平均溫度升高,黃河上游和清水河流域地表蒸發(fā)量明顯增大,汛期降水呈減少趨勢,河道徑流有變小趨勢[22-23]。綜合人類活動和氣候變化因素分析,確定1986年為黃河年最大洪峰序列變異點,1991年為清水河年最大洪峰序列變異點。

    3.2 干支流洪峰時間序列聯(lián)合分布

    圖3 黃河和清水河洪峰發(fā)生時間概率密度擬合

    基于特定EM算法,進行帶懲罰的極大似然估計進行von Mises分布參數(shù)估計,結(jié)果見表3。根據(jù)已估計得到參數(shù),確定兩河時間序列TH和TQ的概率密度函數(shù),并進行擬合,如圖3所示,可以看出,TH和TQ的密度函數(shù)分別選擇4峰和3峰混合分布均能很好擬合經(jīng)驗直方圖。由密度函數(shù)確定分布函數(shù),計算TH和TQ的理論概率,并與經(jīng)驗概率進行擬合,如圖4所示,可以看出,TH和TQ的理論概率均能很好地擬合經(jīng)驗概率。

    圖4 黃河和清水河洪峰發(fā)生時間理論頻率和經(jīng)驗頻率擬合

    河 流α1μ1k1α2μ2k2α3μ3k3α4μ4k4黃 河0.190.551.40.52.853.390.264.627.050.055.8826.04清水河0.041.5170.742.993.90.224.117.91

    由相關(guān)性指標法計算得連接參數(shù)θ=1.034 3,聯(lián)合分布函數(shù)為

    式中:FTH(t)、FTQ(t)分別為TH和TQ的邊緣分布。K-S檢驗法計算聯(lián)合分布統(tǒng)計量D=0.068 4,小于臨界值0.173,通過顯著性檢驗,因此選擇的聯(lián)合分布是合理的。

    3.3 干支流洪峰流量序列聯(lián)合分布

    根據(jù)洪峰序列非一致性檢驗結(jié)果,黃河和清水河洪峰量級序列各存在一個較強變異點,基于可靠性和易于計算原則,對原序列兩部分進行混合分布擬合。由改進模擬退火算法進行參數(shù)估計,結(jié)果見表4。由式(5)確定兩河洪峰量級分布函數(shù),計算理論頻率,與經(jīng)驗頻率進行擬合,并用K-S檢驗法對擬合情況作檢驗,計算的黃河和清水河統(tǒng)計量D分別為0.076和0.126,均小于顯著性水平0.182,構(gòu)建的混合分布能較好地擬合黃河和清水河洪峰量級序列。

    表4 黃河和清水河洪峰量級混合P-Ⅲ分布參數(shù)

    構(gòu)建干支流洪峰發(fā)生時間和量級的部分嵌套四維Gumbel Copula聯(lián)合分布如下:

    C1(TH,QH)=exp{-[(-lnTH)θ1+(-lnQH)θ1]1/θ1}

    (16)

    C2(TQ,QQ)=exp{-[(-lnTQ)θ2+(-lnQQ)θ2]1/θ2}

    (17)

    C(TH,QH,TQ,QQ;θ1,θ2,θ3)=exp{-[(-lnC1)θ3+

    (-lnC2)θ3]1/θ3}

    (18)

    用偽極大似然法估計聯(lián)合分布參數(shù),連接參數(shù)為θ1=1.246,θ2=1.115,θ3=1.065。首先分別對C1和C2進行擬合檢驗,用式(13)(14)計算黃河和清水河各自洪峰發(fā)生時間和量級的兩變量理論聯(lián)合頻率和經(jīng)驗聯(lián)合頻率并進行擬合;進一步計算C1和C2的K-S統(tǒng)計量D1=0.051和D2=0.054,均小于臨界值0.173,通過顯著性檢驗。對C進行擬合檢驗,其K-S統(tǒng)計量D=0.083,小于臨界值0.178,通過顯著性檢驗。這表明上面基于Gumbel Copula構(gòu)建的部分嵌套聯(lián)合分布函數(shù)能很好地擬合黃河和清水河洪峰發(fā)生時間和量級。

    表6 兩河不同頻率洪峰遭遇概率峰值和對應(yīng)的遭遇流量

    3.4 洪水遭遇風(fēng)險分析

    式(10)計算的黃河和清水河洪峰在每一天遭遇的概率如圖5(a)所示。由圖5(a)可以看出,6月15日之前,二者遭遇的概率幾乎為零;7月15日至8月15日為二者主要遭遇時段,也是黃河和清水河洪峰出現(xiàn)的主要時段,在7月31日二者遭遇的概率的達到最大,為3.36×10-4;8月18日至8月31日遭遇概率逐漸增大,并在8月31日出現(xiàn)第二峰值1.15×10-4,主要是由于黃河汛期較長,除7月中下旬至8月上旬的主汛期外,還存在8月末至9月初的次主汛期,而清水河汛期相對滯后于黃河,使得二者在8月25日至9月5日出現(xiàn)第二次遭遇概率高峰。

    式(11)計算的同頻率洪峰在每天遭遇的概率如圖5(b)所示,為便于對比,已將概率值作負對數(shù)處理。從圖5(b)可以看出,兩河同頻率洪峰在7月28日和9月2日附近出現(xiàn)兩個遭遇高峰,表5列出6種重現(xiàn)期同頻率洪峰遭遇的兩個概率峰值;汛初和汛末洪峰遭遇概率相對很小,與黃河和清水河年最大洪峰極少出現(xiàn)在汛初和汛末的實際情況相符。表6為兩河不同頻率洪峰遭遇概率峰值及對應(yīng)的洪峰流量(即青銅峽水庫入庫洪峰流量)。根據(jù)可能發(fā)生遭遇的聯(lián)合洪峰流量,結(jié)合青銅峽水庫調(diào)度規(guī)則,提前預(yù)留足夠防洪庫容,以保證能進行有效攔洪削峰。從表6可以看出,黃河和清水河十年一遇洪峰遭遇的概率為650.91×10-8,未達到萬分之一,且在58 a的實際監(jiān)測中兩河十年一遇洪峰也未發(fā)生遭遇,因此,對下游防洪要求較低的工程,可以忽略二者的遭遇風(fēng)險。

    圖5 黃河和清水河洪峰遭遇概率

    表5 同頻率洪峰遭遇概率峰值

    4 結(jié) 論

    a. 由于氣候、環(huán)境變化和人類活動影響,黃河洪水序列和清水河洪水序列發(fā)生顯著變化的時間不一致,經(jīng)三階段診斷,表明兩河洪水序列發(fā)生顯著變異的年份分別為1986年和1991年。

    b. 不考慮洪峰量級情況下,黃河與清水河洪水發(fā)生時間出現(xiàn)遭遇的概率分別在7月31日和8月31日達到峰值,概率值分別為3.36×10-4和1.15×10-4。

    c. 黃河與清水河同頻率洪峰遭遇概率分別在7月28日和9月2日達到兩個峰值;兩河不同頻率洪峰遭遇概率在7月26日至8月2日和8月31至9月5日兩個時段內(nèi)達到峰值。洪峰遭遇概率均小于萬分之一,因此在實際工程中可以忽略遭遇風(fēng)險。

    [ 1 ] SALVADORI G,MICHELE C D.On the use of Copulas in hydrology: theory and practice[J].Journal of Hydrologic Engineering,2007,12(4):369-380.

    [ 2 ] YAN B W.Coincidence probability of precipitation for the middle route of South-To-North Water Transfer Project in China[J].Journal of Hydrology,2013,499:19-26.

    [ 3 ] GENEST C,FAVRE A C.Meta-elliptical Copulas and their use in frequency analysis of multivariate hydrological data[J].Water Resource Research,2007,43:1-12.

    [ 4 ] CHEN L,SINGH V P.Flood coincidence risk analysis using multivariate Copula functions[J].Journal of Hydrologic Engineering,2012,17(6):742-755.

    [ 5 ] ZHANG Q,GU X H.Flood frequency analysis with consideration of hydrological alterations: changing properties,causes and implications[J].Journal of Hydrology,2014,519:803-813.

    [ 6 ] ALILA Y,MTIRAOUI A.Implications of heterogeneous flood-frequency distributions on traditional stream-discharge prediction techniques[J].Hydrological Processes,2002,16(5):1065-1084.

    [ 7 ] PROHASKA S,ILIC A.Multiple-coincidence of flood waves on the main river and its tributaries[J].IOP Conference Series: Earth and Environmental Science,2008,4:1-8.

    [ 8 ] 謝平,陳廣才,夏軍.變化環(huán)境下基于趨勢分析的水資源評價方法[J].水力發(fā)電學(xué)報,2009,28(2):14-19.(XIE Ping,CHEN Guangcai,XIA Jun.The assessment method of water resources based on trend analysis of changing environments[J].Journal of Hydroelectric Engineering,2009,28(2):14-19.(in Chinese))

    [ 9 ] 謝平,雷紅富,陳廣才,等.水文變異診斷系統(tǒng)[J].水力發(fā)電學(xué)報,2010,29(1).85-90.(XIE Ping,LEI Fuhong,CHEN Guangcai,et al.Hydrological alteration diagnosis system[J].Journal of Hydroelectric Engineering,2010,29(1):85-90.(in Chinese))

    [10] PETTIT A N.A nonparametric approach to the change point problem[J].Applied Statistic 1979,28:126-135.

    [11] PATRICK D B,SHLOMO S S.Increasing physicians awareness of impact of statistics on research outcome: comparative power of thet-test and Wilcoxon Rank-Sum test in small samples applied research[J].Journal of Clinical Epidemiology,1999,52(3):229-253.

    [12] BROWN M B.FORSYTHE A B.The small sample behavior of some statistics which test the equality of several means[J].Techno Metrics,1974,16(1):129-132.[13] 謝平,陳廣才,雷紅富,等.變化環(huán)境下地表水資源評價方法[M].北京:科學(xué)出版社,2009.

    [14] 李瑞雪,張明軍,金爽,等.烏魯木齊河流域氣候變化的區(qū)域差異特征及突變分析[J].干旱區(qū)地理,2010,33(2):243-250.(LI Ruixue,ZHANG Mingjun,JIN Shuang,et al.Regional difference and catastrophe of climate change over Urumqi River Basin[J].Arid Land Geography,2010,33(2):243-250.(in Chinese))

    [15] BRUNSDO C,CORCORAN J.Using circular statistics to analyze time patterns in crime incidence[J].Computers,Environment and Urban Systems,2006,30:300-319.

    [16] FU Y J.Modified likelihood ratio test for homogeneity in a mixture of von Mises distributions[J].Journal of Statistical Planning and Inference,2008,138:667-681.

    [17] 陳家驊,李鵬飛,譚鮮明.混合von Mises模型的參數(shù)估計[J].系統(tǒng)科學(xué)與數(shù)學(xué),2007,27(1):59-67.(CHEN Jiahua,LI Pengfei,TAN Xianming.Inference for von Mises Mixture in mean direction and concentration parameters[J].Journal of System Science and Mathematical Science,2007,27(1): 59-67.(in Chinese))

    [18] CALDERARA S,PRATI A.Mixtures of von Mises distributions for people trajectory shape analysis[J].Transactions on Circuits and Systems for Video Technology,2011,21(4):457-471.

    [19] 成靜清,宋松柏.基于混合分布非一致性年徑流序列頻率參數(shù)的計算[J].西北農(nóng)林科技大學(xué)(自然科學(xué)版),2010,38(2);229-234.(CHENG Jingqing,SONG Songbai.Calculation of hsdrological frequency parameters of inconsistent annual runoff series based on mixed distribution[J].Journal of Northwest Agriculture and Forestry University(Science and Technology),2010,38(2):229-234.(in Chinese))

    [20] 閆寶偉,郭生練,陳璐,等.長江和清江洪水遭遇風(fēng)險分析[J].水利學(xué)報,2010,41(5):553-559.(YAN Baowei,GUO Shenglian,CHEN Lu,et al.Flood encountering risk analysis for the Yangtze River and Qing River[J].Journal of Hydraulic Engineering,2010,41(5):553-559.(in Chinese))

    [21] LIEBSCHER E.Construction of asymmetric multivariate copulas[J].Journal of Multivariate Analysis,2008,99:2234-2250.

    [22] 徐宗學(xué),張楠.黃河流域近50年降水變化趨勢分析[J].地理研究,2006,25(1):27-34.(XU Zongxue,ZHANG Nan.Long term trend of precipitation in the Yellow River Basin during the past 50 years[J].Geographical Research,2006,25(1):27-34.(in Chinese))

    [23] 時興合,秦寧生,汪青春,等.黃河上游徑流變化特征及其影響因素初步分析[J].中國沙漠,2007,27(4):690-697.(SHI Hexing,QIN Ningsheng,WANG Qingchun,et al.Analysis on runoff variation characteristics and influencing factors in the upper Yellow River[J].Journal of Desert Research,2007,27(4):690-697.(in Chinese))

    [24] GENEST C,GHOUDI K,RIVEST L P.A semi-parametric estimation procedure of dependence parameters in multivariate families of distributions[J].Biometrika,1995,85(3):543-552.

    [25] 朱顥東,鐘勇.一種改進的模擬退火算法[J].計算機技術(shù)與發(fā)展,2009,19(6):32-35.(ZHU Haodong,ZHONG Yong.A kind of renewed simulated annealing algorithm[J].Computer Technology and Development,2009,19(6):32-35.(in Chinese))

    Coincidence risk analysis for non-stationary flood peak of Yellow River and its tributaries in Ningxia Hui Autonomous Region

    //LI Ziyuan, FENG Ping, YUAN Ximin

    (StateKeyLaboratoryofHydraulicEngineeringSimulationandSafety,TianjinUniversity,Tianjin300072,China)

    The coincidence risk of the non-stationary flood peak between the Yellow River and Qingshui River was analyzed. The variability of flood peak time series was detected using the comprehensive diagnosis method. The mixed Von Mises distribution and mixed P-Ⅲ distribution were used to fit the flood occurrence dates and peak magnitudes, respectively. The asymmetric partially nested Archimedean Copula function was selected to establish the joint distribution of flood peaks of the two rivers, and the occurrence dates and corresponding magnitudes between the two rivers were studied. Diagrams of coincidence risk of floods with different combinations of frequencies of the two rivers were provided. The results indicate that the flood time series in both rivers had changed; the coincidence risk of the flood peak occurrence dates of the two rivers peaked on July 31 and August 31, with values of 3.36×10-4and 1.15×10-4, respectively; and the peak coincidence risks of flood peaks with difference frequencies were obtained in two periods, from July 26 to August 2 and from August 31 to September 5.

    flood peak coincidence; non-stationary; mixed distribution; Copula function; risk analysis; Yellow River; Qingshui River

    10.3880/j.issn.1006-7647.2016.06.010

    國家自然科學(xué)基金(51179117,51209157)

    李子遠(1991—),男,碩士研究生,主要從事水文及水資源研究。E-mail: 362919252@qq.com

    TV122+.5

    A

    1006-7647(2016)06-0051-07

    2015-10-10 編輯:鄭孝宇)

    猜你喜歡
    兩河清水河洪峰
    清水河邊
    飛天(2022年5期)2022-05-18 08:11:45
    古代兩河流域文明運河功能探析
    清水河生態(tài)清潔小流域
    陸西地區(qū)清水河組一段儲層特征及差異性分析
    淡定!
    解禁洪峰
    “一江兩河”區(qū)域青稞氮肥推薦指標體系研究
    西藏科技(2016年9期)2016-09-26 12:21:35
    兩河流域早王朝時期作為地理概念的“蘇美爾”
    清水河上游流域可收集雨水資源量估算與檢驗
    地火(2014年4期)2014-03-01 01:55:30
    国产黄片美女视频| 日韩精品青青久久久久久| 岛国在线观看网站| 亚洲国产高清在线一区二区三 | 午夜免费成人在线视频| 久久精品夜夜夜夜夜久久蜜豆 | 亚洲国产毛片av蜜桃av| 日韩成人在线观看一区二区三区| 久久人人精品亚洲av| 午夜福利一区二区在线看| 88av欧美| 久久中文字幕一级| 丝袜人妻中文字幕| 久久精品亚洲精品国产色婷小说| 成人精品一区二区免费| 激情在线观看视频在线高清| 黄网站色视频无遮挡免费观看| 午夜视频精品福利| 啦啦啦免费观看视频1| 岛国视频午夜一区免费看| 女性生殖器流出的白浆| 久久人妻福利社区极品人妻图片| cao死你这个sao货| 亚洲五月天丁香| 亚洲一卡2卡3卡4卡5卡精品中文| 少妇的丰满在线观看| 麻豆国产av国片精品| 最近在线观看免费完整版| 国产成人欧美在线观看| 精品国产乱子伦一区二区三区| 一级a爱视频在线免费观看| www.自偷自拍.com| 啦啦啦免费观看视频1| 亚洲av电影在线进入| 国产国语露脸激情在线看| 久久 成人 亚洲| 国产精品综合久久久久久久免费| 国产精品日韩av在线免费观看| 亚洲一卡2卡3卡4卡5卡精品中文| 91在线观看av| 国产一区在线观看成人免费| 午夜福利在线观看吧| svipshipincom国产片| 亚洲av五月六月丁香网| 在线免费观看的www视频| www.www免费av| 国产av一区在线观看免费| 亚洲av五月六月丁香网| 亚洲精品色激情综合| 日韩大码丰满熟妇| 99精品在免费线老司机午夜| ponron亚洲| 日韩三级视频一区二区三区| 欧美日韩一级在线毛片| 人人妻人人看人人澡| 好男人电影高清在线观看| 国产亚洲欧美在线一区二区| 在线av久久热| 少妇熟女aⅴ在线视频| 亚洲国产精品999在线| 叶爱在线成人免费视频播放| 色av中文字幕| 亚洲五月色婷婷综合| 精品国产美女av久久久久小说| 18美女黄网站色大片免费观看| 亚洲第一电影网av| 侵犯人妻中文字幕一二三四区| 国产av一区二区精品久久| 久久天躁狠狠躁夜夜2o2o| 成人18禁在线播放| 国产精品一区二区精品视频观看| 99精品欧美一区二区三区四区| 1024手机看黄色片| 一级毛片精品| 99国产极品粉嫩在线观看| 国产精品精品国产色婷婷| 两人在一起打扑克的视频| 国产成+人综合+亚洲专区| 日韩成人在线观看一区二区三区| 欧美zozozo另类| 熟女少妇亚洲综合色aaa.| 亚洲国产精品合色在线| 少妇被粗大的猛进出69影院| 99riav亚洲国产免费| 两个人看的免费小视频| 日韩精品青青久久久久久| av在线天堂中文字幕| 首页视频小说图片口味搜索| 国产视频一区二区在线看| 精品乱码久久久久久99久播| 久久久精品国产亚洲av高清涩受| 欧美日韩瑟瑟在线播放| 99热只有精品国产| 欧美色欧美亚洲另类二区| 欧美 亚洲 国产 日韩一| 丁香欧美五月| 亚洲成av人片免费观看| 中文资源天堂在线| av在线播放免费不卡| 亚洲五月天丁香| 欧美久久黑人一区二区| 欧美一区二区精品小视频在线| 午夜福利免费观看在线| 校园春色视频在线观看| 日本成人三级电影网站| 久久久久久国产a免费观看| 免费在线观看亚洲国产| 在线观看66精品国产| 免费在线观看完整版高清| 欧美人与性动交α欧美精品济南到| 国产精品野战在线观看| 国产精品香港三级国产av潘金莲| 桃红色精品国产亚洲av| 国产亚洲精品一区二区www| 日韩成人在线观看一区二区三区| 午夜福利在线观看吧| 久久热在线av| 久久人人精品亚洲av| 人成视频在线观看免费观看| 免费观看精品视频网站| 国产成人系列免费观看| 又大又爽又粗| 丁香六月欧美| avwww免费| 十八禁人妻一区二区| 脱女人内裤的视频| 欧美在线黄色| 一级作爱视频免费观看| 天天一区二区日本电影三级| 久久婷婷人人爽人人干人人爱| a在线观看视频网站| 中文字幕久久专区| 精华霜和精华液先用哪个| 亚洲五月天丁香| 亚洲人成网站高清观看| 亚洲自拍偷在线| 最近最新免费中文字幕在线| 久久久久久久久免费视频了| 无遮挡黄片免费观看| 国产精品久久久久久人妻精品电影| 黄色丝袜av网址大全| 国产高清视频在线播放一区| 亚洲av第一区精品v没综合| 正在播放国产对白刺激| 99国产极品粉嫩在线观看| 国产成人精品久久二区二区91| 亚洲欧洲精品一区二区精品久久久| 亚洲精品国产区一区二| 18禁美女被吸乳视频| 精品久久久久久久久久久久久 | 欧美成人性av电影在线观看| 人妻丰满熟妇av一区二区三区| 久久久久亚洲av毛片大全| 亚洲一码二码三码区别大吗| 国产精品野战在线观看| 男女做爰动态图高潮gif福利片| 黄片大片在线免费观看| 热re99久久国产66热| 丁香欧美五月| 欧美日韩一级在线毛片| av免费在线观看网站| 亚洲国产欧洲综合997久久, | 日韩欧美一区视频在线观看| 在线观看日韩欧美| 国产精品98久久久久久宅男小说| 亚洲人成77777在线视频| 欧美日韩亚洲国产一区二区在线观看| 欧美精品啪啪一区二区三区| 哪里可以看免费的av片| 熟女少妇亚洲综合色aaa.| 国产一级毛片七仙女欲春2 | 亚洲一区二区三区色噜噜| 午夜福利免费观看在线| 制服人妻中文乱码| 午夜福利在线在线| 亚洲五月婷婷丁香| 岛国在线观看网站| 久久久水蜜桃国产精品网| 久久精品国产99精品国产亚洲性色| 搡老熟女国产l中国老女人| 少妇被粗大的猛进出69影院| 国产亚洲欧美98| 少妇 在线观看| 琪琪午夜伦伦电影理论片6080| 人妻久久中文字幕网| 精品一区二区三区视频在线观看免费| 欧美日韩中文字幕国产精品一区二区三区| 正在播放国产对白刺激| 一级a爱片免费观看的视频| 午夜免费激情av| 色综合婷婷激情| 久久99热这里只有精品18| 在线观看免费午夜福利视频| 又紧又爽又黄一区二区| 亚洲av第一区精品v没综合| 九色国产91popny在线| 免费看a级黄色片| 丝袜人妻中文字幕| www日本在线高清视频| 国产成年人精品一区二区| 亚洲免费av在线视频| 亚洲五月色婷婷综合| 久久久精品欧美日韩精品| 国产高清videossex| 日韩欧美国产一区二区入口| 成人手机av| 香蕉久久夜色| 欧美激情极品国产一区二区三区| 天堂影院成人在线观看| 国产在线观看jvid| 亚洲成人久久爱视频| 一本大道久久a久久精品| 欧美黄色片欧美黄色片| 日韩欧美在线二视频| 好男人电影高清在线观看| 亚洲精华国产精华精| 成人精品一区二区免费| 亚洲国产欧洲综合997久久, | 午夜免费观看网址| 十分钟在线观看高清视频www| 中文亚洲av片在线观看爽| 国产精品自产拍在线观看55亚洲| 欧美最黄视频在线播放免费| 他把我摸到了高潮在线观看| 成年免费大片在线观看| 国内少妇人妻偷人精品xxx网站 | 免费女性裸体啪啪无遮挡网站| 国产在线精品亚洲第一网站| 黄色女人牲交| 人妻久久中文字幕网| 国产爱豆传媒在线观看 | 高清在线国产一区| 精品国产一区二区三区四区第35| 午夜激情福利司机影院| 麻豆久久精品国产亚洲av| 日日干狠狠操夜夜爽| 精品国产乱子伦一区二区三区| 午夜福利欧美成人| 人人妻人人澡人人看| 变态另类成人亚洲欧美熟女| 狠狠狠狠99中文字幕| 亚洲精品av麻豆狂野| 久久精品国产亚洲av香蕉五月| 99久久国产精品久久久| 亚洲成人久久性| 99国产精品一区二区蜜桃av| av电影中文网址| 天堂动漫精品| 亚洲欧洲精品一区二区精品久久久| 国产野战对白在线观看| 亚洲成av人片免费观看| 一边摸一边做爽爽视频免费| 亚洲国产精品成人综合色| 国产1区2区3区精品| 少妇被粗大的猛进出69影院| 国产亚洲欧美精品永久| 波多野结衣av一区二区av| 99精品久久久久人妻精品| 亚洲精品一卡2卡三卡4卡5卡| 天天添夜夜摸| 成人精品一区二区免费| 日韩三级视频一区二区三区| 黄色丝袜av网址大全| 日韩成人在线观看一区二区三区| 日韩一卡2卡3卡4卡2021年| 亚洲三区欧美一区| 最近最新免费中文字幕在线| 手机成人av网站| 嫩草影院精品99| 黄片播放在线免费| 大型黄色视频在线免费观看| av有码第一页| 久久久国产成人免费| 女人爽到高潮嗷嗷叫在线视频| 一区福利在线观看| 久久精品国产99精品国产亚洲性色| 午夜福利高清视频| 一级作爱视频免费观看| 999精品在线视频| 久久亚洲真实| 色老头精品视频在线观看| a级毛片a级免费在线| 夜夜躁狠狠躁天天躁| 一区二区日韩欧美中文字幕| 国产亚洲欧美精品永久| 人人妻人人澡欧美一区二区| 热99re8久久精品国产| 久久精品91蜜桃| 日韩av在线大香蕉| 久久精品国产综合久久久| 一二三四社区在线视频社区8| 精品久久久久久成人av| 国产真人三级小视频在线观看| 亚洲人成伊人成综合网2020| 此物有八面人人有两片| 国产极品粉嫩免费观看在线| 一进一出抽搐gif免费好疼| 视频在线观看一区二区三区| 国产精品一区二区精品视频观看| 啪啪无遮挡十八禁网站| 婷婷亚洲欧美| 大型av网站在线播放| 亚洲一区二区三区不卡视频| 久久久国产精品麻豆| 夜夜躁狠狠躁天天躁| 亚洲精品久久成人aⅴ小说| 国产成人欧美在线观看| 精品久久久久久,| 人人妻人人澡人人看| 丝袜在线中文字幕| 久久久国产成人免费| 亚洲熟女毛片儿| 国产精品 国内视频| 在线观看舔阴道视频| 黄色视频不卡| 免费一级毛片在线播放高清视频| 国产三级黄色录像| 久久 成人 亚洲| 亚洲国产欧洲综合997久久, | videosex国产| 精品高清国产在线一区| 久久久久久久午夜电影| 看片在线看免费视频| 一二三四社区在线视频社区8| 天堂影院成人在线观看| 亚洲成人国产一区在线观看| 国产又黄又爽又无遮挡在线| 午夜免费观看网址| 久久国产精品人妻蜜桃| 亚洲第一青青草原| 黑人欧美特级aaaaaa片| 欧美午夜高清在线| 国产av一区在线观看免费| 日韩精品免费视频一区二区三区| 成熟少妇高潮喷水视频| 国产精品香港三级国产av潘金莲| 法律面前人人平等表现在哪些方面| 久久中文字幕人妻熟女| 欧美乱码精品一区二区三区| 久久中文字幕人妻熟女| 欧美乱色亚洲激情| 男女之事视频高清在线观看| 女性被躁到高潮视频| 亚洲精品粉嫩美女一区| 啦啦啦韩国在线观看视频| 国产黄色小视频在线观看| 男男h啪啪无遮挡| 女生性感内裤真人,穿戴方法视频| 91老司机精品| 变态另类丝袜制服| 精品久久久久久久久久免费视频| 亚洲精品美女久久久久99蜜臀| 中文字幕最新亚洲高清| 精品乱码久久久久久99久播| 高潮久久久久久久久久久不卡| 久久中文看片网| 国产又黄又爽又无遮挡在线| 亚洲成av片中文字幕在线观看| 黑人欧美特级aaaaaa片| 亚洲精品美女久久av网站| 日韩大尺度精品在线看网址| 亚洲精品粉嫩美女一区| 色综合婷婷激情| 90打野战视频偷拍视频| 国产亚洲精品一区二区www| 极品教师在线免费播放| 不卡av一区二区三区| 制服丝袜大香蕉在线| 久久精品亚洲精品国产色婷小说| 两个人视频免费观看高清| 国产免费男女视频| 国产成人精品久久二区二区91| 少妇熟女aⅴ在线视频| av视频在线观看入口| 欧美精品啪啪一区二区三区| 一级片免费观看大全| 成年版毛片免费区| 18禁裸乳无遮挡免费网站照片 | 69av精品久久久久久| 欧美激情久久久久久爽电影| 18美女黄网站色大片免费观看| 一级黄色大片毛片| 香蕉国产在线看| 亚洲精品国产区一区二| 国产国语露脸激情在线看| 日韩精品青青久久久久久| 免费在线观看亚洲国产| 在线视频色国产色| 一卡2卡三卡四卡精品乱码亚洲| 亚洲国产看品久久| 精品久久久久久久久久久久久 | 精品久久久久久成人av| 国产一卡二卡三卡精品| 精品国产乱子伦一区二区三区| 久久久水蜜桃国产精品网| 免费一级毛片在线播放高清视频| 极品教师在线免费播放| 色老头精品视频在线观看| 亚洲国产欧美一区二区综合| 国产精品99久久99久久久不卡| 麻豆av在线久日| 亚洲人成网站在线播放欧美日韩| 香蕉国产在线看| ponron亚洲| 很黄的视频免费| 久久人人精品亚洲av| 99久久综合精品五月天人人| 色综合亚洲欧美另类图片| 久久国产精品人妻蜜桃| 国产一区二区在线av高清观看| 亚洲av中文字字幕乱码综合 | 成年版毛片免费区| 黄片大片在线免费观看| 亚洲全国av大片| 夜夜夜夜夜久久久久| 国产成人系列免费观看| 亚洲第一电影网av| 一二三四社区在线视频社区8| 免费一级毛片在线播放高清视频| 国内毛片毛片毛片毛片毛片| 老司机在亚洲福利影院| 欧美乱妇无乱码| 国产爱豆传媒在线观看 | 久久精品亚洲精品国产色婷小说| 国产精品久久久久久人妻精品电影| 在线观看66精品国产| 成人手机av| 999久久久国产精品视频| 99久久精品国产亚洲精品| 精华霜和精华液先用哪个| 日本三级黄在线观看| 岛国视频午夜一区免费看| 可以在线观看毛片的网站| 亚洲专区中文字幕在线| 久久久久国内视频| 999久久久精品免费观看国产| 99精品久久久久人妻精品| 欧美又色又爽又黄视频| 国产亚洲欧美精品永久| 亚洲,欧美精品.| 亚洲 欧美一区二区三区| 嫩草影院精品99| 亚洲最大成人中文| 男人的好看免费观看在线视频 | 亚洲成av片中文字幕在线观看| 久久 成人 亚洲| 日韩欧美免费精品| 婷婷精品国产亚洲av| 黄片大片在线免费观看| 97碰自拍视频| 亚洲av电影不卡..在线观看| 美女高潮喷水抽搐中文字幕| 精品久久久久久久人妻蜜臀av| 少妇裸体淫交视频免费看高清 | 老熟妇仑乱视频hdxx| av片东京热男人的天堂| 一区福利在线观看| 啦啦啦观看免费观看视频高清| 90打野战视频偷拍视频| 曰老女人黄片| 一本一本综合久久| 人人妻,人人澡人人爽秒播| 成人三级做爰电影| 亚洲,欧美精品.| 一级毛片高清免费大全| 日本a在线网址| 欧美成人午夜精品| 非洲黑人性xxxx精品又粗又长| www日本黄色视频网| 久久午夜综合久久蜜桃| videosex国产| 手机成人av网站| 成年人黄色毛片网站| 99国产极品粉嫩在线观看| 欧美日韩亚洲国产一区二区在线观看| 99精品欧美一区二区三区四区| 亚洲av片天天在线观看| 母亲3免费完整高清在线观看| 欧美日韩黄片免| 日韩一卡2卡3卡4卡2021年| 一a级毛片在线观看| 亚洲人成网站高清观看| 久久中文字幕人妻熟女| 在线观看免费视频日本深夜| 国产成人一区二区三区免费视频网站| 男人舔女人下体高潮全视频| 精品久久久久久久久久免费视频| 叶爱在线成人免费视频播放| 女同久久另类99精品国产91| 非洲黑人性xxxx精品又粗又长| 久久狼人影院| 啦啦啦观看免费观看视频高清| 久久中文字幕一级| 侵犯人妻中文字幕一二三四区| 可以在线观看的亚洲视频| 日韩成人在线观看一区二区三区| 久久热在线av| 免费看美女性在线毛片视频| 午夜两性在线视频| av片东京热男人的天堂| 国产欧美日韩一区二区三| 一级毛片精品| 色播亚洲综合网| 99久久久亚洲精品蜜臀av| 久久草成人影院| 深夜精品福利| 亚洲自偷自拍图片 自拍| 这个男人来自地球电影免费观看| 97碰自拍视频| www.熟女人妻精品国产| 露出奶头的视频| 精品福利观看| 亚洲成人久久爱视频| 国产不卡一卡二| 高潮久久久久久久久久久不卡| 免费观看人在逋| 国产亚洲av高清不卡| 精品无人区乱码1区二区| 两个人免费观看高清视频| 脱女人内裤的视频| 久久久久久久精品吃奶| 午夜福利高清视频| 欧洲精品卡2卡3卡4卡5卡区| 色综合站精品国产| 国产精品免费视频内射| 很黄的视频免费| avwww免费| 一级毛片女人18水好多| 国产成人精品久久二区二区91| 丝袜在线中文字幕| 亚洲av电影不卡..在线观看| 欧美激情极品国产一区二区三区| 亚洲av五月六月丁香网| 国产成人av教育| 久久精品aⅴ一区二区三区四区| 激情在线观看视频在线高清| 在线观看免费午夜福利视频| 亚洲片人在线观看| 国产成人精品久久二区二区免费| 视频在线观看一区二区三区| 亚洲真实伦在线观看| 成人18禁高潮啪啪吃奶动态图| АⅤ资源中文在线天堂| 成人国产综合亚洲| 人人妻人人澡人人看| 黄色女人牲交| 在线观看舔阴道视频| 午夜免费激情av| 亚洲欧美日韩无卡精品| 怎么达到女性高潮| 亚洲av熟女| 18禁观看日本| 午夜成年电影在线免费观看| 国产精品爽爽va在线观看网站 | 欧美zozozo另类| 午夜福利高清视频| 丰满人妻熟妇乱又伦精品不卡| 久久久水蜜桃国产精品网| 成年版毛片免费区| 美女 人体艺术 gogo| 国产欧美日韩一区二区三| 美国免费a级毛片| tocl精华| 午夜福利欧美成人| 精品高清国产在线一区| 国产成人影院久久av| 色在线成人网| 看免费av毛片| av超薄肉色丝袜交足视频| 婷婷精品国产亚洲av在线| xxx96com| 大型av网站在线播放| 嫩草影视91久久| 成人18禁高潮啪啪吃奶动态图| 国产黄色小视频在线观看| 国产高清视频在线播放一区| svipshipincom国产片| 黑丝袜美女国产一区| 国产午夜福利久久久久久| 99久久综合精品五月天人人| 美女高潮到喷水免费观看| 777久久人妻少妇嫩草av网站| aaaaa片日本免费| 日本精品一区二区三区蜜桃| 成人国产综合亚洲| 日韩一卡2卡3卡4卡2021年| 国产成年人精品一区二区| 美女国产高潮福利片在线看| 午夜激情福利司机影院| 日韩欧美三级三区| 欧美日韩瑟瑟在线播放| 中文字幕精品亚洲无线码一区 | 日本在线视频免费播放| 久久中文字幕人妻熟女| 99久久久亚洲精品蜜臀av| 国产精品 欧美亚洲| 久久中文字幕人妻熟女| 成人av一区二区三区在线看| svipshipincom国产片| 日韩欧美免费精品| 精品电影一区二区在线| 国产aⅴ精品一区二区三区波| 日日干狠狠操夜夜爽| 男人舔奶头视频| 黄色a级毛片大全视频| 99久久无色码亚洲精品果冻| 2021天堂中文幕一二区在线观 | 精品免费久久久久久久清纯| 大型黄色视频在线免费观看| 久久热在线av| 啪啪无遮挡十八禁网站| 亚洲自拍偷在线| 99国产精品99久久久久| 给我免费播放毛片高清在线观看| 日韩精品免费视频一区二区三区|