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

    基于Copula函數(shù)的松花江流域水文干旱頻率分析

    2021-03-27 07:55:30廖顯薇高峰魏婷宋小燕宋松柏
    灌溉排水學(xué)報 2021年3期
    關(guān)鍵詞:松花江流域佳木斯歷時

    廖顯薇,高峰,魏婷,宋小燕*,宋松柏

    基于Copula函數(shù)的松花江流域水文干旱頻率分析

    廖顯薇1,高峰2,魏婷1,宋小燕1*,宋松柏1

    (1.西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院,陜西 楊凌 712100;2.內(nèi)蒙古自治區(qū)水利科學(xué)研究院,呼和浩特 010051)

    【】基于Copula函數(shù)對松花江流域進行水文干旱頻率分析。利用可變閾值法識別日尺度水文干旱事件的歷時和烈度,并進行融合處理;選擇6種常用的分布函數(shù)擬合干旱特征變量,以Kolmogorov-Smirnov 檢驗法優(yōu)選單變量邊緣分布;從3種Copula函數(shù)中優(yōu)選函數(shù)模型,建立了干旱歷時和烈度二維聯(lián)合分布,并計算聯(lián)合分布概率和重現(xiàn)期,再利用蒙特卡洛方法,討論基于Copula干旱研究的不確定性。融合標準取c=0.1和c=5,流域干旱歷時均值為81~105 d;對數(shù)正態(tài)分布是大賚站、扶余站和哈爾濱站干旱歷時的最優(yōu)邊緣分布,而各個站點干旱烈度最優(yōu)分布各異;Frank Copula為流域干旱二維聯(lián)合分布的最優(yōu)模型,流域大多數(shù)干旱的重現(xiàn)期不超過20a;在同現(xiàn)重現(xiàn)期為20 a的條件下,大賚、扶余、哈爾濱和佳木斯4個水文站點的最大可能設(shè)計值表現(xiàn)出較大的不確定性。Copula函數(shù)能夠較好地擬合松花江流域水文干旱特征變量間的聯(lián)合分布,但需考慮其不確定性影響。

    水文干旱;Copula函數(shù);干旱事件融合;松花江流域

    0 引言

    【研究意義】干旱是一種由水分的收支或供求不平衡而形成的水分短缺現(xiàn)象[1]。隨著全球人口增長、氣候變暖,各地區(qū)對水分的需求量增加,導(dǎo)致干旱頻發(fā)[2]。東北地區(qū)1960—2000年因水旱災(zāi)害糧食損失1 498億kg,直接經(jīng)濟損失達1 585億元[3]。東北地區(qū)是重要的糧食生產(chǎn)基地,對中國的糧食安全起著保障作用。松花江作為黑龍江的支流,在東北地區(qū)承擔(dān)著重要的社會和經(jīng)濟責(zé)任。因此,研究松花江流域的干旱問題,對流域的農(nóng)業(yè)生產(chǎn)和生態(tài)管理等具有重要意義。

    【研究進展】美國氣象學(xué)會將干旱分為氣象干旱、水文干旱、農(nóng)業(yè)干旱及社會經(jīng)濟干旱4個類型[4]?,F(xiàn)階段對松花江流域的干旱研究多集中于氣象干旱,而水文干旱研究不足[5]。水文干旱是由河川徑流或地下水資源不足導(dǎo)致的結(jié)果,相較于氣象和農(nóng)業(yè)干旱,水文干旱的發(fā)生更加復(fù)雜,干旱程度更深、后果更嚴重[6]。構(gòu)建多變量聯(lián)合分布及計算重現(xiàn)期是目前水文干旱研究的重點和難點。Copula函數(shù)可克服傳統(tǒng)多變量分布函數(shù)要求同一類邊緣分布函數(shù)的限制[7],單參數(shù)Achimedean Copula函數(shù)族更因構(gòu)造簡單,相對容易求解,被廣泛應(yīng)用于水文領(lǐng)域[8-10]?!厩腥朦c】干旱特征序列樣本長度較短,代表性不高,基于樣本序列開展的聯(lián)合分布研究存在顯著的不確定性,這對流域水資源合理管理提出了巨大挑戰(zhàn)[11]。Copula函數(shù)雖已在各個流域的水文干旱頻率分析中取得良好的研究成果,但目前關(guān)于Copula函數(shù)的不確定性研究的報道較少?!緮M解決的關(guān)鍵問題】本文選取松花江流域大賚、扶余、哈爾濱和佳木斯4個水文站點1960—2015年的日平均流量,采用可變閾值法識別松花江流域日尺度干旱,通過敏感性測試,確定符合流域水文特性的干旱事件融合標準,再選取常用的伽馬、威布爾、對數(shù)正態(tài)等6種概率分布函數(shù)擬合優(yōu)選干旱特征變量的邊緣分布,在此基礎(chǔ)上利用Copula函數(shù)構(gòu)建二維干旱變量的聯(lián)合分布函數(shù)并計算重現(xiàn)期,對計算結(jié)果進行不確定性探討。

    1 材料與方法

    1.1 研究區(qū)概況

    松花江流域位于中國東北部,有南北兩源,北源為嫩江,南源為第二松花江,流域面積為55.68萬km2,年徑流量為762億m3。流域地處北溫帶季風(fēng)氣候區(qū),大陸性氣候特點非常明顯,冬季寒冷漫長,夏季炎熱多雨,春季干燥多風(fēng),秋季很短,年內(nèi)溫差較大,多年平均氣溫在3~5 ℃之間。流域多年平均降水量一般在500 mm左右,汛期6—9月的降水量占全年60%~80%,冬季12—次年2月的降水量僅為全年的5%左右。由于松花江流域主要以大氣降水補給為主、融雪補給為輔,因此徑流量的年內(nèi)分配也具有明顯的季節(jié)特征。

    本文選取研究區(qū)內(nèi)大賚站(嫩江匯入松花江干流控制站)、扶余站(第二松花江匯入松花江干流控制站)、哈爾濱站(松花江干流中游控制站)和佳木斯站(松花江干流下游控制站)4個水文站點1960—2015年日平均流量。

    1.2 研究方法

    1.2.1 基于可變閾值法識別干旱歷時和烈度

    以4個水文站點的日平均流量數(shù)據(jù)為基礎(chǔ),使用基于30 d滑動窗口流量—持續(xù)時間曲線第80個百分點的可變閾值提取干旱事件的歷時和烈度[12]。這意味著,1年中的每一天都有1個不同的閾值水平,這種方法已被用于季節(jié)性顯著的集水區(qū)[13]。當(dāng)流量低于閾值時,記為干旱,干旱的時長為干旱歷時;低于閾值的總?cè)彼繛楦珊盗叶取T谔崛〉?個歷時較長且烈度較大的相鄰水文干旱事件之間,存在某一個日徑流量超過設(shè)定閾值的過程,該過程歷時較短,烈度較小,通常認為其前后2個干旱事件是相互關(guān)聯(lián)的,因此,需要建立1個相鄰干旱事件的融合標準。以時間間隔t和間隔時間內(nèi)超閾徑流量i與前一個干旱事件所缺徑流量的比值i為判別標準[14],當(dāng)相鄰干旱事件的tp都低于臨界值c和c時,則有:

    式中:{d,s}和{d1,s1}為2個相鄰干旱事件;p和p為相鄰干旱事件融合后的干旱歷時和烈度。

    通過敏感性測試[15],c取不同值(0、0.05、0.1、0.15、0.2、0.25);對于日徑流過程,c一般取5 d左右(c取2、3、4、5、6 d),對原干旱特征序列(c=0)進行融合,將融合后的歷時與烈度均值與原序列均值做比較,得出均值比隨c和c的變化曲線。當(dāng)多條曲線斜率同時變緩時確定c的值,曲線間落差較大時確定c的值,則可得出融合標準。以哈爾濱站為例,其敏感性測試見圖1和圖2。由圖1和圖2可以看出,當(dāng)c=0.1時,5條曲線斜率明顯變小,說明干旱歷時和烈度的均值變化趨于平緩;當(dāng)c=6時,干旱均值產(chǎn)生大的跳躍。大賚站、扶余站和佳木斯站均有此情況,因此,本文取c=0.1和c=5兩個融合標準。

    圖1 哈爾濱站干旱歷時敏感性曲線

    圖2 哈爾濱站干旱烈度敏感性曲線

    1.2.2 干旱歷時和干旱烈度的邊緣分布

    基于Copula函數(shù)的多變量頻率分析主要分為單變量頻率分析和多變量Copula頻率分析2部分[16],確定干旱特征變量的邊緣分布是進行Copula函數(shù)選擇和參數(shù)計算的關(guān)鍵步驟。受氣候變化和高強度人類活動影響,水文事件的概率分布函數(shù)未知,可通過選取大量候選概率分布函數(shù)和擬合優(yōu)度檢驗,篩選其中的最優(yōu)者。本文采用常用的伽馬、威布爾、對數(shù)正態(tài)等6種概率分布函數(shù)來擬合干旱歷時和烈度,以Kolmogorov-Smirnov檢驗(K-S)進行擬合優(yōu)選。本文中邊緣分布及聯(lián)合分布參數(shù)均采用最大似然法估計。

    1.2.3 基于Copula函數(shù)的聯(lián)合分布

    采用Copula函數(shù)構(gòu)建干旱特征變量聯(lián)合分布之前,需確定干旱特征變量間的相依性。本文采用Kendall、Spearman秩相關(guān)系數(shù)和Pearson古典相關(guān)系數(shù)來描述干旱歷時和烈度間的相依性。

    Copula函數(shù)是定義域為[0,1]的均勻分布的多維聯(lián)合分布函數(shù),它可以將多個隨機變量的邊緣分布連接起來得到它們的聯(lián)合分布。

    本文選取常用的3種單參數(shù)Archimedean Copula函數(shù)(Gumbel Copula、Clayton Copula和Frank Copula)構(gòu)建干旱特征變量的聯(lián)合分布,其函數(shù)和參數(shù)取值范圍見表1。

    表1 Copulas函數(shù)及參數(shù)取值范圍

    采用準則和均方根誤差檢驗Copula函數(shù)的擬合優(yōu)度,檢驗標準值越小,說明函數(shù)的擬合程度越好。

    1.2.4 干旱歷時和干旱烈度的二維聯(lián)合分布與重現(xiàn)期

    1)二維聯(lián)合概率分布

    假定干旱歷時的邊緣分布函數(shù)為D(),干旱烈度的邊緣分布函數(shù)為S(),干旱歷時和烈度的聯(lián)合分布為:

    (,)=(≤,≤)=[D(),S()]=(,)。 (2)

    干旱歷時和烈度聯(lián)合超越概率為:

    給定≥時,的條件概率:

    2)重現(xiàn)期

    聯(lián)合重現(xiàn)期:

    同現(xiàn)重現(xiàn)期:

    式中:為系列長度;為時段內(nèi)干旱發(fā)生次數(shù)。

    1.2.5 基于蒙特卡洛的聯(lián)合概率計算不確定性評估

    蒙特卡洛方法作為一種常見的分析水文數(shù)值模擬不確定性方法,具有實用性廣、方法簡便的特點,且可將參數(shù)的不確定性轉(zhuǎn)化為模擬結(jié)果的不確定性[17]。本文利用蒙特卡洛方法,對樣本進行多重抽樣,計算Copula聯(lián)合分布下給定重現(xiàn)期的最大可能設(shè)計值[18]的置信區(qū)間,置信區(qū)間能夠用來描述水文設(shè)計值的不確定性,提供相對可靠的變動范圍,可評估樣本不確定性對干旱聯(lián)合分布的影響,詳細步驟見文獻[19]。

    2 結(jié)果與分析

    2.1 干旱特征變量分析

    根據(jù)可變閾值法提取干旱特征變量并進行統(tǒng)計可知,1960—2015年,松花江流域水文干旱事件年均發(fā)生3.27~4.16次,水文干旱次數(shù)發(fā)生最多的為中游哈爾濱站,其次為第二松花江扶余站和下游佳木斯站,最少為嫩江大賚站,各水文站點干旱歷時和烈度的統(tǒng)計特征(最大值、上四分位數(shù)3、均值、中位數(shù)、下四分位數(shù)1、最小值)見表2。

    對干旱歷時進行分析可知,大賚站、扶余站、哈爾濱站和佳木斯站歷時大于30 d的干旱次數(shù)分別占干旱總次數(shù)的89%、73%、84%和86%;其中歷時大于1 a的干旱分別發(fā)生2、4、1和3次;歷時最長的干旱分別發(fā)生在2012、2014、1965年和1984年;4個水文站的干旱歷時主要集中在30~90 d。

    表2 水文干旱歷時和烈度統(tǒng)計特征

    大賚站年平均徑流量為199.67億m3,扶余站年平均徑流量為145.22億m3,哈爾濱站年平均徑流量為396.42億m3,佳木斯站年平均徑流量為614.51億m3。4個站點干旱烈度(缺水量)與年平均徑流之比大于10%的干旱發(fā)生次數(shù)分別為8、4、2和3次,其中烈度最大的干旱分別發(fā)生在1998、2005、1998年和1984年,與年平均徑流量之比分別為:50%、37%、34%和35%。各個站點的年缺水量與年徑流量如圖3所示。由圖3可知,4個站點的年徑流量均呈下降趨勢,年缺水量隨徑流量的增減而波動,說明采用可變閾值法識別的干旱烈度與徑流量關(guān)聯(lián)性很強,除佳木斯站有略微下降趨勢外,其余3站均有略微上升趨勢,大賚站的干旱情況最嚴重。在干旱事件中,大賚站=28 d,=23.21億m3、扶余站=25 d,=35.63億m3、哈爾濱站=26 d,=36.82億m3和佳木斯站=20 d,=30.25億m3等干旱事件,雖然干旱歷時都小于1個月,但干旱烈度較大,這些“短強干旱”事件大多發(fā)生在夏季,其次是春季和秋季,冬季最少。

    對干旱歷時與烈度序列進行相關(guān)性分析,結(jié)果見表3。

    從表3可以得出,干旱烈度和歷時間存在較強的相關(guān)性。因此,可以采用Copula函數(shù)建立兩變量的聯(lián)合分布模型。

    2.2 邊緣分布函數(shù)的擬合

    采用伽馬(GAMA)、威布爾(WBL)、對數(shù)正態(tài)(LOGN)、指數(shù)(EXP)、正態(tài)(NOR)和廣義極值(GEV)等6種概率分布函數(shù)來擬合干旱歷時和烈度,以K-S方法進行擬合優(yōu)選,結(jié)果見表4。

    表4 水文干旱特征值邊緣分布模型擬合優(yōu)選

    從表4可以得出,LOGN分布為大賚站、扶余站和哈爾濱站干旱歷時的最優(yōu)邊緣分布模型,3個站點干旱烈度的最優(yōu)邊緣分布模型分別為GAMA模型、WBL模型和GEV模型。佳木斯站干旱歷時和烈度的最優(yōu)邊緣分布模型為GEV模型和EXP模型。

    2.3 Copula函數(shù)優(yōu)選

    利用3種Copula函數(shù)構(gòu)建干旱歷時和烈度的二維聯(lián)合分布,以和對模型進行擬合優(yōu)度評價,結(jié)果見表5。由表5可知,F(xiàn)rank Copula對大賚站、扶余站、佳木斯站和哈爾濱站的干旱歷時與烈度之間的擬合效果最好,優(yōu)選Frank Copula建立4個站點的二維聯(lián)合概率及重現(xiàn)期計算。

    表5 Copulas函數(shù)擬合優(yōu)度評價

    注 標*數(shù)據(jù)對應(yīng)的和值最小。

    2.4 干旱歷時和烈度的聯(lián)合概率

    4個水文站點的干旱歷時和烈度超越聯(lián)合概率以及條件概率分布規(guī)律相類似,聯(lián)合概率值隨干旱歷時和烈度值的減小而呈現(xiàn)增大趨勢,表明易發(fā)生歷時較短且干旱烈度較小的干旱。由圖4可以得到,不同歷時(≥)下的干旱烈度(≤)所發(fā)生的概率。因篇幅限制,圖4以佳木斯站為例進行規(guī)律闡述。當(dāng)干旱烈度不大于某一值時,其條件概率隨歷時的增大而減小,例如當(dāng)≤50億m3時,(≥31)>(≥61)>(≥91)。

    圖4 佳木斯站干旱歷時與烈度的條件概率

    2.5 重現(xiàn)期

    利用Frank Copula構(gòu)建聯(lián)合分布,再根據(jù)式(5)和式(6)計算4個站點的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期,以各個站點干旱歷時與烈度的上四分位數(shù)統(tǒng)計特征值為界,流域大多數(shù)干旱的重現(xiàn)期不超過20 a。大賚站和扶余站的聯(lián)合重現(xiàn)期比哈爾濱站和佳木斯站小,說明嫩江和第二松花江在遭遇長干旱歷時或強干旱烈度的干旱事件的概率要大于松花江干流。而同現(xiàn)重現(xiàn)期則相反,特別是哈爾濱站的同現(xiàn)重現(xiàn)期有內(nèi)凹的趨勢,容易發(fā)生歷時短且烈度小的干旱。說明流域面積越大,對徑流的調(diào)節(jié)能力越強,對“或”事件的抵御能力越強[20]。

    本文重點分析同現(xiàn)重現(xiàn)期為20 a所對應(yīng)的最大可能設(shè)計值,采用蒙特卡洛方法,對干旱特征序列進行多重抽樣,獲得給定條件下的區(qū)間估計,從而討論其不確定性,聯(lián)合分布的置信度取95%。分析得出,同現(xiàn)重現(xiàn)期隨著干旱歷時和烈度的增大而增大。當(dāng)同現(xiàn)重現(xiàn)期為20a時,所有水文站點的聯(lián)合設(shè)計值二元置信區(qū)間面積跨度都較大,哈爾濱站和佳木斯站的置信區(qū)間可從5a跨越到50a。將干旱歷時與烈度聯(lián)合設(shè)計值二元置信與其單變量邊緣分布結(jié)合來看,4個水文站點的干旱歷時與烈度邊緣分布曲線都出現(xiàn)了因水文數(shù)據(jù)不滿足正態(tài)分布假設(shè),具有的尖峰、重尾特征,這影響對極端干旱事件的估計,導(dǎo)致干旱歷時與烈度聯(lián)合設(shè)計值二元置信區(qū)間面積大,增加聯(lián)合設(shè)計值的不確定性。在4個水文站點中,扶余站和大賚站的聯(lián)合設(shè)計值二元置信區(qū)間面積跨度沒有超過50a重現(xiàn)期,相對松花江干流的哈爾濱站和佳木斯站不確定性較小,原因可能是位于第二松花江干流的控制性工程豐滿水庫和白山水庫與位于嫩江干流的控制性工程尼爾基水庫對兩個區(qū)域的水文干旱有調(diào)節(jié)作用,可以減少聯(lián)合設(shè)計值的不確定性。

    3 討論

    本研究采用可變閾值法,能夠適應(yīng)徑流季節(jié)性變化,通過敏感性試驗,選取的融合標準與涂新軍等[21]在東江流域的研究一致,研究方法可靠。松花江流域的水文干旱歷時均值3個月,水文干旱有加重趨勢,嫩江大賚站的干旱情況最嚴重,與韓冬梅等[22]的研究結(jié)果相符。除此之外,值得注意的是,研究結(jié)果中一些“短強干旱”事件具有明顯的季節(jié)性特征,大多發(fā)生在夏季,其次是春季和秋季。這應(yīng)與松花江流域以降水補給為主、融雪補給為輔的不均勻來水情況有關(guān)。夏秋季遇強降雨時會形成洪澇,春季受融雪補給的作用,各河流會形成程度不等的春汛,當(dāng)短期內(nèi)發(fā)生洪水或豐枯季節(jié)交替時,日流量變異大,導(dǎo)致與閾值的差積增大。不僅如此,徑流年內(nèi)分配的不均勻性在區(qū)域上也存在差異。相較于第二松花江和松花江干流,嫩江流域的水文站點年內(nèi)徑流分配最不均勻[23],這可能是導(dǎo)致大賚站干旱情況嚴重的原因之一。

    在獲取干旱特征的基礎(chǔ)上,本文優(yōu)選Frank Copula構(gòu)建聯(lián)合分布模型并計算松花江流域水文干旱重現(xiàn)期,重點討論同現(xiàn)重現(xiàn)期為20a的條件下,最大可能設(shè)計值的二元置信區(qū)間,表現(xiàn)出很強的不確定性,大賚站和扶余站的不確定性較松花江干流小。水文設(shè)計值不確定的影響因素主要集中在3個方面[24]:①水文事件本身的不確定性;②線型選擇的不確定性;③模型參數(shù)估計方法的不確定性。松花江干流有松嫩平原和三江平原,居民數(shù)量多,活動范圍廣,人類活動通過改變大氣循環(huán)和通過蓄、引、提、調(diào)水工程等改變河流蓄存狀態(tài)和水力聯(lián)系[25],加大了水文事件的不確定性;干旱特征邊緣分布線型的選取缺乏物理根據(jù),僅僅是一種數(shù)值擬合,因此線型選擇的偏差,也存在不確定性[26];分布函數(shù)確定情況下,參數(shù)估計方法不同,其水文設(shè)計值置信區(qū)間估計方法不同。李航等[27]研究6種參數(shù)估計方法發(fā)現(xiàn),優(yōu)化適線法和線性矩法在水文頻率分析中較極大似然法取得的擬合效果更為理想。綜上所述,聯(lián)合分布不確定性是多因素所致的客觀現(xiàn)象,今后的研究可通過對初始條件的確定、模型的優(yōu)化等途徑,使輸出結(jié)果更加精確[28]。

    已有很多研究將Copula函數(shù)用于各個流域的水文干旱頻率分析計算中,在前人研究的基礎(chǔ)上,本文討論了基于Copula函數(shù)聯(lián)合分布的不確定性,對流域管理決策的制定有非常重要的意義。本文也存在一些不足,如干旱特征除歷時和烈度外,還有峰值和面積等特征,二維聯(lián)合分布仍不能滿足實際需要;其次,在干旱年際對比上,可變閾值法表現(xiàn)出一定的不足,它展現(xiàn)不出豐枯年份的差異,且未將干旱烈度標準化,無法直接體現(xiàn)干旱嚴重的空間異質(zhì)性。另外,本文僅對水文干旱頻率的不確定性進行了研究,如何解決不確定性產(chǎn)生的影響等問題,需在今后的研究中進行更深入的探討。目前Copula函數(shù)已有上百種,如何將其成熟的用于干旱研究仍是未來一段時間內(nèi)的研究難點。

    4 結(jié)論

    1)本文運用可變閾值法對松花江流域4個水文站的日徑流量進行干旱識別,取c=5和c=0.1為融合標準,得到松花江流域水文干旱歷時均值為81~105 d,烈度均值為8.52億~33.68億m3,且存在歷時小于1個月,烈度較大的干旱事件。

    2)對數(shù)正態(tài)、伽馬、威布爾和廣義極值等分布能較好地擬合干旱歷時和烈度的邊緣分布;Frank Copula對大賚站、扶余站、哈爾濱站和佳木斯站的干旱歷時與烈度之間的擬合效果最好。

    3)超越聯(lián)合概率值隨干旱歷時和烈度的減小而呈增大趨勢,而重現(xiàn)期則相反。松花江流域所發(fā)生的干旱大多數(shù)重現(xiàn)期不超過20a。在同現(xiàn)重現(xiàn)期為20a的條件下,4個站點的最大可能設(shè)計值置信區(qū)間較大,最大跨越5~50a,不利于流域工程設(shè)計值的確定。

    [1] 洪興駿, 郭生練, 王樂, 等. 基于最大熵原理的水文干旱指標計算方法研究[J]. 南水北調(diào)與水利科技, 2018, 16(2): 93-99.

    HONG Xingjun, GUO Shenglian, WANG Le, et al. Hydrological drought index calculation using Principle of Maximum Entropy[J]. South-to-North Water Transfers and Water Science & Technology, 2018, 16(2): 93-99.

    [2] 趙福年, 王潤元, 王鶯, 等. 干旱過程?時空尺度及干旱指數(shù)構(gòu)建機制的探討[J]. 災(zāi)害學(xué), 2018, 33(4): 32-39.

    ZHAO Funian, WANG Runyuan, WANG Ying, et al. Discussion of process, temporal and spatial for drought and establishment of drought indexs[J]. Journal of Catastrophology, 2018, 33(4): 32-39.

    [3] 水利部松遼水利委員會. 東北區(qū)水旱災(zāi)害集[M]. 長春: 吉林出版社, 2003: 25-81.

    Songliao Water Resources Commission, Ministry of Water Resources. Flood and drought disasters in Northeast China[M]. Changchun: Jilin Press, 2003: 25-81.

    [4] HEIM J, RICHARD R. A review of twentieth-century drought Indices used in the united states[J]. Bulletin of the American Meteorological Society, 2002, 83(8): 1 149-1 165.

    [5] 吳燕鋒, 章光新. 松花江區(qū)氣象水文干旱演變特征[J]. 地理科學(xué), 2018, 38(10): 1 731-1 739.

    WU Yanfeng, ZHANG Guangxin. Spatio-temporal patterns of meteorological and hydrological drought in the Songhua River Area from 1961 to 2010[J]. Scientia Geographica Sinica, 2018, 38(10): 1 731-1 739.

    [6] 董前進, 謝平. 水文干旱研究進展[J]. 水文, 2014, 34(4): 1-7.

    DONG Qianjin, XIE Ping. Advances in hydrological research[J]. Hydrology, 2014, 34(4): 1-7.

    [7] 宋松柏, 王小軍. 基于Copula函數(shù)的水文隨機變量和概率分布計算[J]. 水利學(xué)報, 2018, 49(6): 687-693.

    SONG Songbai, WANG Xiaojun. Probability distribution calculation of the sum of hydrological random variables based on Copula function approach[J]. Journal of Hydraulic Engineering, 2018, 49(6): 687-693.

    [8] 袁超, 宋松柏. 基于Copula函數(shù)的水文干旱聯(lián)合概率分布研究[J]. 水利與建筑工程學(xué)報, 2016, 14(6): 50-53.

    YUAN Chao, SONG Songbai. The bivariate hydrologic drought distribution based on copula function[J]. Journal of Water Resources and Architectural Engineering, 2016, 14(6): 50-53.

    [9] 張迎, 黃生志, 黃強, 等. 基于Copula函數(shù)的新型綜合干旱指數(shù)構(gòu)建與應(yīng)用[J]. 水利學(xué)報, 2018, 49(6): 703-714.

    ZHANG Ying, HUANG Shengzhi, HUANG Qiang, et al. Construction and application of a new comprehensive drought index based on Copula function[J]. Journal of Hydraulic Engineering, 2018, 49(6): 703-714.

    [10] 顧世祥, 趙眾, 陳晶, 等. 基于高維Copula函數(shù)的逐日潛在蒸散量及氣象干旱預(yù)測[J]. 農(nóng)業(yè)工程學(xué)報, 2020, 36(9): 143-151.

    GU Shixiang, ZHAO Zhong, CHEN Jing, et al. Daily reference evapotranspiration and meteorological drought forecast using high-dimensional Copula joint distribution model[J]. Chinese Society of Agricultural Engineering, 2020, 36(9): 143-151.

    [11] ZHANG Q, XIAO M Z, SINGH V P. Uncertainty evaluation of copula analysis of hydrological droughts in the East River basin, China[J]. Global and Planetary Change, 2015, 129(6): 1-9.

    [12] BEYENE B, VAN LOON A F, VAN LANEN H A J, et al. Investigation of variable threshold level approaches for hydrological drought identification[J]. Journal of Hydrology, 2015, 512(11): 421-434.

    [13] VAN LOON A F, LAAHA G. Hydrological drought severity explained by climate and catchment characteristics[J]. Journal of Hydrology, 2015, 526(12): 3-14.

    [14] MADSEN H, ROSBJERG D. On the modelling of extreme droughts[J]. Water Resources Systems, 1995, 231(6): 377-385.

    [15] LENA M T, HENRIK M, BENTE C. On the definition and modelling of streamflow drought duration and deficit volume[J]. Hydrological Sciences Journal, 1997, 42(1): 15-33.

    [16] 徐翔宇, 許凱, 楊大文, 等. 多變量干旱事件識別與頻率計算方法[J]. 水科學(xué)進展, 2019, 30(3): 373-381.

    XU Xiangyu, XU Kai, YANG Dawen, et al. Drought identification and drought frequency analysis based on multiple variables[J]. Advances in Water Science, 2019, 30(3): 373-381.

    [17] 苗添升, 盧文喜, 歐陽琦, 等. 地下水?dāng)?shù)值模擬的不確定性分析在水質(zhì)預(yù)測中的應(yīng)用[J]. 水電能源科學(xué), 2016, 34(8): 20-23, 44.

    MIAO Tiansheng, LU Wenxi, OUYANG Qi, et al. Application of Uncertainty Analysis of Groundwater Numerical Simulation in Water Quality Prediction[J]. Water Resources and Power, 2016, 34(8): 20-23, 44.

    [18] SALVADORI G, DE MICHELE C, DURANTE F. On the return period and design in a multivariate framework[J]. Hydrology and Earth System Sciences, 2011, 15(11): 3 293-3 305.

    [19] 馬川惠, 黃強, 郭愛軍. 涇河流域水沙聯(lián)合分布特征分析及其不確定性評估[J]. 水利學(xué)報, 2019, 50(2): 273-282.

    MA Chuanhui, HUANG Qiang, GUO Aijun. Characteristic analysis and uncertainty assessment of joint distribution of flow and sand in Jinghe River basin[J]. Journal of Hydraulic Engineering, 2019, 50(2): 273-282.

    [20] 姚蕊, 孫鵬, 張強, 等. 基于Copula淮河流域水文干旱頻率分析及影響研究[J]. 北京師范大學(xué)學(xué)報(自然科學(xué)版), 2019, 55(6): 755-763.

    YAO Rui, SUN Peng, ZHANG Qiang, et al. Hydrological drought frequency based on Copulas in Huai River Basin[J]. Journal of Beijing Normal University Natural Sciences, 2019, 55(6): 755-763.

    [21] 涂新軍, 陳曉宏, 趙勇, 等. 變化環(huán)境下東江流域水文干旱特征及缺水響應(yīng)[J]. 水科學(xué)進展, 2016, 27(6): 810-821.

    TU Xinjun, CHEN Xiaohong, ZHAO Yong, et al. Responses of hydrological drought properties and water shortage under changing environments in Dongjiang River basin[J]. Advances in Water Science, 2016, 27(6): 810-821.

    [22] 韓冬梅, 楊貴羽, 嚴登華, 等. 近50年東北地區(qū)旱澇時空特征分析[J]. 水電能源科學(xué), 2014, 32(6): 5-8.

    HAN Dongmei, YANG Guiyu, YAN Denghua, et al. Spatial-temporal feature analysis of drought and flood in Northeast China in recent 50 years[J]. Water Water Resources and Power, 2014, 32(6): 5-8.

    [23] 李峰平. 變化環(huán)境下松花江流域水文與水資源響應(yīng)研究[D]. 長春: 中國科學(xué)院研究生院(東北地理與農(nóng)業(yè)生態(tài)研究所), 2015.

    LI Fengping. Study on water cycle and water resources in Songhua River Basin under the Changing Environment[D]. Changchun: Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, 2015.

    [24] 雷冠軍, 王文川, 殷峻暹, 等. P-Ⅲ型曲線參數(shù)估計方法研究綜述[J].人民黃河, 2017, 39(10): 1-7.

    LEI Guanjun, WANG Wenchuan, YIN Junxian, et al. Review on study of parameter estimation method on P-Ⅲcurve[J]. Yellow River, 2017, 39(10): 1-7.

    [25] 王文, 王靖淑, 陶奕源, 等. 人類活動對水文干旱形成與發(fā)展的影響研究進展[J]. 水文, 2020, 40(3): 1-8.

    WANG Wen, WANG Jingshu, TAO Yiyuan, et al. Review of human impacts on formation and development of hydrologic drought[J]. Journal of China Hydrology, 2020, 40(3): 1-8.

    [26] 劉攀, 郭生練, 田向榮, 等. 基于貝葉斯理論的水文頻率線型選擇與綜合[J]. 武漢大學(xué)學(xué)報(工學(xué)版), 2005, 38(5): 38-42.

    LIU Pan, Guo Shenglian, Tian Xiangrong, et al. Selecting and averaging of flood frequency models based on Bayesian theory[J]. Engineering Journal of Wuhan University, 2005, 38(5): 38-42.

    [27] 李航, 宋松柏, 石繼海. 指數(shù)Gamma分布參數(shù)估計方法對比研究[J]. 水力發(fā)電學(xué)報, 2019, 38(4): 96-107.

    LI Hang, SONG Songbai, SHI Jihai. Comparison of parameter estimation methods for exponential Gamma distribution[J]. Jourmal of Hydroelectric Engineering, 2019, 38(4): 96-107.

    [28] 強安豐, 汪妮, 莫淑紅, 等. 氣候變化對水文水資源影響評價的不確定研究進展[J].水資源研究, 2020, 9(2): 169-178.

    QIANG Anfeng, WANG Ni, MO Shuhong, et al. Review for impact assessment of climate change on hydrology and water resources in uncertainties research[J]. Journal of Water Resources Research, 2020, 9(2): 169-178.

    Using Copula Method to Analyze Drought Frequency in Songhua River Basin

    LIAO Xianwei1, GAO Feng2, WEI Ting1, SONG Xiaoyan1*, SONG Songbai1

    (1. College of Water Resources and Architectural Engineering, Northwest A&F University, Yangling 712100, China 2.Inner Mongolia Water Conservancy Research Institute, Hohhot 010051, China)

    Drought could significantly impact hydrological processes and water resources, and understanding frequency and uncertainty of drought occurrence is hence imperative to water resources planning. The copula method has emerged as an improved multivariate analysis over the univariate analysis for quantitative analysis of drought.【】This paper is to present the results of the frequency and severity of droughts in Songhua river basin calculated using the copula multivariate method.【】The analysis was based on archived data, and duration and severity of the pooled daily drought events were identified using the variable threshold level method. The drought index variables were fitted by six common distribution functions respectively, and they were then evaluated using the Kolmogorov Smirnov test method. Optimal model selected from three Copula functions was used toestablish a two-dimensional joint distribution for the drought index variables, from which we calculated the joint distribution probability as well as the return period.The uncertainty of the drought was analyzed using the Monte Carlo method.【】With the pooling criteria set to bec=0.1 andc=5, the average drought duration was 81~105 days in the basin. The lognormal distribution can adequately describe the marginal distribution of the drought duration in Dalai, Fuyu and Harbin stations in the basin, despite the differences in optimal distribution of the drought severity between them. The Frank copula was the best model for the two-dimensional joint distribution of the drought, and the return periods of the droughts were less than 20 years. When the return period was 20 years, the maximum design values for Dalai, Fuyu, Harbin and Jiamusi hydrological stations were most uncertain.【】The copula multivariate model can adequately describe the joint distribution of the drought index variables in Songhua river basin, and its application should consider drought uncertainty.

    hydrologicaldrought; Copula function; pooling of drought events; Songhua river basin

    P333

    A

    10.13522/j.cnki.ggps.2020330

    1672 - 3317(2021)03 - 0134 - 08

    廖顯薇,高峰, 魏婷, 等. 基于Copula函數(shù)的松花江流域水文干旱頻率分析[J]. 灌溉排水學(xué)報,2021, 40(3): 134-141.

    LIAO Xianwei, GAO Feng, WEI Ting, et al. Using Copula Method to Analyze Drought Frequency in Songhua River Basin [J]. Journal of Irrigation and Drainage, 2021, 40(3): 134-141.

    2020-06-19

    國家自然科學(xué)基金項目(41501022);國家科技基礎(chǔ)資源調(diào)查專項(2017FY100904);中央高?;究蒲袠I(yè)務(wù)項目(2452020167)

    廖顯薇(1998-),女。碩士研究生,主要從事水文干旱研究。E-mail: 1070928072@qq.com

    宋小燕(1983-),女。副教授,博士,主要從事水文循環(huán)方面的研究。E-mail: xiaoyansong@nwsuaf.edu.cn

    責(zé)任編輯:韓 洋

    猜你喜歡
    松花江流域佳木斯歷時
    自然災(zāi)害、農(nóng)戶傳統(tǒng)風(fēng)險管理手段與天氣指數(shù)保險需求
    淺談佳木斯站減速頂?shù)膽?yīng)用、管理及設(shè)備的專業(yè)化維修
    量詞“只”的形成及其歷時演變
    常用詞“怠”“惰”“懶”的歷時演變
    佳木斯防爆電機研究所
    防爆電機(2020年6期)2020-12-14 07:16:52
    佳木斯大學(xué)張莉莉水彩作品選登
    吉林省松花江流域產(chǎn)業(yè)生態(tài)系統(tǒng)適應(yīng)性分析
    對《紅樓夢》中“不好死了”與“……好的”的歷時考察
    古今字“兌”“說”“悅”“敚”歷時考察
    基于一維水動力模型的洪水頂托影響分析
    日韩国内少妇激情av| 久久亚洲真实| 久久久国产精品麻豆| 美女午夜性视频免费| 亚洲精品中文字幕在线视频| 久久国产亚洲av麻豆专区| 怎么达到女性高潮| 欧洲精品卡2卡3卡4卡5卡区| 一本久久中文字幕| 亚洲av成人av| 18禁观看日本| 国产精品亚洲av一区麻豆| 国产一区二区三区综合在线观看| 欧美精品啪啪一区二区三区| 成人手机av| 亚洲精品久久成人aⅴ小说| 成人国产一区最新在线观看| 午夜免费成人在线视频| 搡老岳熟女国产| 久久国产精品男人的天堂亚洲| av中文乱码字幕在线| av中文乱码字幕在线| 午夜福利18| 无遮挡黄片免费观看| 久久人人精品亚洲av| 成人三级黄色视频| 一级a爱片免费观看的视频| 欧美在线黄色| 国产精品av久久久久免费| 99精品久久久久人妻精品| 中文字幕最新亚洲高清| 在线观看66精品国产| 成人三级黄色视频| 亚洲av第一区精品v没综合| 非洲黑人性xxxx精品又粗又长| 两个人看的免费小视频| 黄色成人免费大全| 久久狼人影院| 18禁美女被吸乳视频| 国产麻豆成人av免费视频| 窝窝影院91人妻| 亚洲电影在线观看av| 在线永久观看黄色视频| 狠狠狠狠99中文字幕| 好男人在线观看高清免费视频 | 国产成人av教育| 免费不卡黄色视频| 啦啦啦韩国在线观看视频| 精品久久久精品久久久| 一级黄色大片毛片| 少妇粗大呻吟视频| 国产在线观看jvid| 两个人免费观看高清视频| 久久精品aⅴ一区二区三区四区| 午夜福利,免费看| or卡值多少钱| 国产av又大| 欧美成人免费av一区二区三区| 一区二区三区国产精品乱码| 精品少妇一区二区三区视频日本电影| 午夜福利18| 色综合站精品国产| 男人的好看免费观看在线视频 | 国产精品香港三级国产av潘金莲| 精品国产美女av久久久久小说| 黄色视频,在线免费观看| 免费在线观看亚洲国产| 免费不卡黄色视频| 国产激情欧美一区二区| 国产精品亚洲一级av第二区| 午夜影院日韩av| 成人三级黄色视频| 国产精品久久久久久精品电影 | av有码第一页| 丝袜美腿诱惑在线| 一个人观看的视频www高清免费观看 | 男女下面插进去视频免费观看| 老司机福利观看| 一边摸一边抽搐一进一出视频| 啦啦啦免费观看视频1| av视频在线观看入口| 自拍欧美九色日韩亚洲蝌蚪91| 日韩精品中文字幕看吧| 欧美午夜高清在线| 在线观看免费日韩欧美大片| 免费久久久久久久精品成人欧美视频| 国产成人一区二区三区免费视频网站| 高潮久久久久久久久久久不卡| 黄片大片在线免费观看| 视频区欧美日本亚洲| avwww免费| 国产91精品成人一区二区三区| 国产私拍福利视频在线观看| 久久婷婷成人综合色麻豆| 欧美日本中文国产一区发布| 老司机深夜福利视频在线观看| 国产国语露脸激情在线看| 国产午夜精品久久久久久| 好男人电影高清在线观看| 9191精品国产免费久久| 丝袜美足系列| 丝袜美足系列| 亚洲一区中文字幕在线| 不卡av一区二区三区| 91九色精品人成在线观看| 欧美成人午夜精品| 99久久久亚洲精品蜜臀av| 亚洲欧美精品综合一区二区三区| 日韩国内少妇激情av| 少妇粗大呻吟视频| 午夜福利一区二区在线看| 一个人免费在线观看的高清视频| 精品久久久久久久人妻蜜臀av | 国产亚洲精品久久久久久毛片| 日本 av在线| 亚洲激情在线av| 国产欧美日韩精品亚洲av| 99精品欧美一区二区三区四区| 欧美日韩中文字幕国产精品一区二区三区 | 国产免费av片在线观看野外av| 男人的好看免费观看在线视频 | 亚洲天堂国产精品一区在线| 12—13女人毛片做爰片一| videosex国产| 精品久久蜜臀av无| 日韩中文字幕欧美一区二区| 成在线人永久免费视频| 国产av一区在线观看免费| 妹子高潮喷水视频| 久久中文字幕人妻熟女| 亚洲国产中文字幕在线视频| 欧美一区二区精品小视频在线| 亚洲av电影在线进入| 亚洲欧美精品综合久久99| 亚洲专区国产一区二区| 国产欧美日韩精品亚洲av| 极品人妻少妇av视频| 757午夜福利合集在线观看| 国产精品久久久久久精品电影 | 国产1区2区3区精品| 国产野战对白在线观看| 又大又爽又粗| 欧美乱码精品一区二区三区| 黄片小视频在线播放| 欧美日韩福利视频一区二区| 国产伦一二天堂av在线观看| 国产私拍福利视频在线观看| 97超级碰碰碰精品色视频在线观看| 欧美精品啪啪一区二区三区| 国产黄a三级三级三级人| 一进一出抽搐动态| 国产精品1区2区在线观看.| av欧美777| 国产精品av久久久久免费| 久久香蕉激情| 多毛熟女@视频| 丝袜美足系列| 国产主播在线观看一区二区| 青草久久国产| 欧美精品啪啪一区二区三区| 美女高潮到喷水免费观看| 首页视频小说图片口味搜索| 欧美日韩黄片免| 亚洲精品粉嫩美女一区| 青草久久国产| 免费在线观看黄色视频的| 欧美一级毛片孕妇| 一级毛片精品| 又黄又粗又硬又大视频| 成人特级黄色片久久久久久久| 大型黄色视频在线免费观看| av在线播放免费不卡| 久久精品aⅴ一区二区三区四区| 亚洲在线自拍视频| 免费久久久久久久精品成人欧美视频| 免费观看精品视频网站| 精品久久久精品久久久| 亚洲精品在线美女| 日韩 欧美 亚洲 中文字幕| 久久久国产精品麻豆| 欧美成人性av电影在线观看| 亚洲av第一区精品v没综合| 亚洲成a人片在线一区二区| 丝袜美足系列| 日韩免费av在线播放| 亚洲久久久国产精品| 两个人视频免费观看高清| 嫁个100分男人电影在线观看| www.精华液| 精品少妇一区二区三区视频日本电影| 在线观看一区二区三区| 欧美激情极品国产一区二区三区| 欧美人与性动交α欧美精品济南到| 国产亚洲精品一区二区www| 少妇被粗大的猛进出69影院| 在线av久久热| 亚洲一区中文字幕在线| 国产高清视频在线播放一区| 亚洲第一欧美日韩一区二区三区| 亚洲天堂国产精品一区在线| 757午夜福利合集在线观看| 村上凉子中文字幕在线| 如日韩欧美国产精品一区二区三区| 亚洲情色 制服丝袜| 国产乱人伦免费视频| АⅤ资源中文在线天堂| 校园春色视频在线观看| 伊人久久大香线蕉亚洲五| 欧美日韩亚洲综合一区二区三区_| 精品熟女少妇八av免费久了| 亚洲欧美精品综合一区二区三区| 亚洲av成人av| 精品国产一区二区久久| 国产99白浆流出| 久久人妻熟女aⅴ| 黄片播放在线免费| 亚洲欧美日韩无卡精品| 无人区码免费观看不卡| 欧美日韩黄片免| 亚洲七黄色美女视频| 亚洲情色 制服丝袜| 欧美日本亚洲视频在线播放| 久9热在线精品视频| 99国产精品一区二区三区| 亚洲av成人av| 成人18禁在线播放| 巨乳人妻的诱惑在线观看| 一级作爱视频免费观看| 黑人巨大精品欧美一区二区mp4| 国产一区二区三区在线臀色熟女| 淫妇啪啪啪对白视频| 97人妻精品一区二区三区麻豆 | cao死你这个sao货| 亚洲一区二区三区不卡视频| 十八禁网站免费在线| 校园春色视频在线观看| 亚洲少妇的诱惑av| 中文字幕人妻熟女乱码| 欧美日韩精品网址| 久久精品人人爽人人爽视色| 50天的宝宝边吃奶边哭怎么回事| 国产精品久久久人人做人人爽| av片东京热男人的天堂| 大型黄色视频在线免费观看| 后天国语完整版免费观看| 日韩精品中文字幕看吧| 免费高清在线观看日韩| 久久 成人 亚洲| 两个人免费观看高清视频| 欧美日韩福利视频一区二区| 亚洲av成人一区二区三| 女警被强在线播放| 制服人妻中文乱码| 国产精品免费视频内射| 精品久久久久久,| 女性生殖器流出的白浆| 一个人观看的视频www高清免费观看 | 国产精品乱码一区二三区的特点 | 丁香六月欧美| 两个人免费观看高清视频| 女人精品久久久久毛片| 乱人伦中国视频| 久久九九热精品免费| 成人国产综合亚洲| 51午夜福利影视在线观看| 中文字幕av电影在线播放| 国产欧美日韩一区二区三| 可以免费在线观看a视频的电影网站| 国产一级毛片七仙女欲春2 | 国内精品久久久久精免费| 日本 av在线| 极品人妻少妇av视频| 国产一级毛片七仙女欲春2 | 91老司机精品| 国产一卡二卡三卡精品| 在线播放国产精品三级| 精品一品国产午夜福利视频| 人妻久久中文字幕网| 婷婷精品国产亚洲av在线| 欧美色视频一区免费| 国产精品99久久99久久久不卡| 黑丝袜美女国产一区| 怎么达到女性高潮| 天堂动漫精品| 最新美女视频免费是黄的| 国产又色又爽无遮挡免费看| 中出人妻视频一区二区| 午夜影院日韩av| 久99久视频精品免费| 午夜两性在线视频| 欧美成人一区二区免费高清观看 | 给我免费播放毛片高清在线观看| 国产精品98久久久久久宅男小说| 最近最新中文字幕大全免费视频| 搡老熟女国产l中国老女人| 中文字幕人妻丝袜一区二区| 正在播放国产对白刺激| 成人三级黄色视频| 狂野欧美激情性xxxx| 此物有八面人人有两片| 日本一区二区免费在线视频| 丝袜在线中文字幕| 国产欧美日韩一区二区精品| 999精品在线视频| 丝袜人妻中文字幕| 黑人巨大精品欧美一区二区mp4| 亚洲一码二码三码区别大吗| 亚洲最大成人中文| 精品免费久久久久久久清纯| 国产免费男女视频| 精品国产乱码久久久久久男人| 免费一级毛片在线播放高清视频 | 国产成人精品久久二区二区免费| 久久香蕉国产精品| 久久 成人 亚洲| 99国产精品免费福利视频| 中文字幕最新亚洲高清| 性色av乱码一区二区三区2| 日韩大尺度精品在线看网址 | 精品一区二区三区av网在线观看| 久久中文看片网| 国产亚洲欧美精品永久| 村上凉子中文字幕在线| 女人精品久久久久毛片| 午夜影院日韩av| 国产成人av教育| 悠悠久久av| 91国产中文字幕| 久久国产精品影院| 日韩欧美免费精品| www.精华液| 一级黄色大片毛片| 国产xxxxx性猛交| 老鸭窝网址在线观看| 人人澡人人妻人| 欧美日韩精品网址| 国产一卡二卡三卡精品| 满18在线观看网站| 校园春色视频在线观看| 天天一区二区日本电影三级 | 日本vs欧美在线观看视频| 99国产精品一区二区三区| 免费少妇av软件| 国产精品久久久久久人妻精品电影| av片东京热男人的天堂| 好看av亚洲va欧美ⅴa在| 国产亚洲精品第一综合不卡| 久久久久久久午夜电影| 精品无人区乱码1区二区| 亚洲国产欧美网| 国产主播在线观看一区二区| av福利片在线| 免费无遮挡裸体视频| 亚洲精品国产精品久久久不卡| 久久精品国产99精品国产亚洲性色 | 超碰成人久久| 日韩 欧美 亚洲 中文字幕| av超薄肉色丝袜交足视频| 午夜免费鲁丝| 99国产极品粉嫩在线观看| 免费搜索国产男女视频| 精品久久久久久久久久免费视频| 涩涩av久久男人的天堂| 18禁美女被吸乳视频| 午夜免费鲁丝| 成人三级黄色视频| 给我免费播放毛片高清在线观看| 黄色a级毛片大全视频| 搡老岳熟女国产| 欧美乱妇无乱码| 成人国产综合亚洲| 亚洲一码二码三码区别大吗| 亚洲va日本ⅴa欧美va伊人久久| 国产av精品麻豆| 一本综合久久免费| 国产成人精品久久二区二区91| 欧美绝顶高潮抽搐喷水| 满18在线观看网站| 一级a爱视频在线免费观看| 麻豆av在线久日| 少妇 在线观看| 日韩大码丰满熟妇| 国产一级毛片七仙女欲春2 | 国产日韩一区二区三区精品不卡| 国产精品一区二区三区四区久久 | 老司机深夜福利视频在线观看| 国产伦人伦偷精品视频| 曰老女人黄片| 又大又爽又粗| 免费在线观看影片大全网站| 女警被强在线播放| 久久人人精品亚洲av| 少妇的丰满在线观看| 亚洲欧美日韩无卡精品| 亚洲男人的天堂狠狠| 婷婷六月久久综合丁香| 最近最新中文字幕大全电影3 | 精品乱码久久久久久99久播| 国产97色在线日韩免费| 久久久久久人人人人人| 99riav亚洲国产免费| 无限看片的www在线观看| 国产午夜精品久久久久久| 九色国产91popny在线| 久久国产精品男人的天堂亚洲| 99久久国产精品久久久| 亚洲精品国产区一区二| 久久久久久久久久久久大奶| 一本大道久久a久久精品| 岛国视频午夜一区免费看| 成人三级做爰电影| 亚洲欧美日韩无卡精品| av片东京热男人的天堂| 亚洲中文字幕一区二区三区有码在线看 | 欧美日韩中文字幕国产精品一区二区三区 | 一级毛片女人18水好多| 国产私拍福利视频在线观看| 大码成人一级视频| 黄网站色视频无遮挡免费观看| 性色av乱码一区二区三区2| 夜夜躁狠狠躁天天躁| 色在线成人网| 午夜精品在线福利| 欧美精品亚洲一区二区| 亚洲精品美女久久久久99蜜臀| 欧美中文综合在线视频| 国产三级在线视频| 青草久久国产| 欧美日韩瑟瑟在线播放| 国产成年人精品一区二区| 很黄的视频免费| 老司机靠b影院| 国产精品久久视频播放| 欧美 亚洲 国产 日韩一| 国产成人精品在线电影| 人妻丰满熟妇av一区二区三区| 精品国产美女av久久久久小说| 国产精品一区二区免费欧美| 午夜福利免费观看在线| 两性午夜刺激爽爽歪歪视频在线观看 | 手机成人av网站| 制服丝袜大香蕉在线| 成人国产综合亚洲| 亚洲五月婷婷丁香| 日韩精品青青久久久久久| 首页视频小说图片口味搜索| 日韩欧美一区二区三区在线观看| 精品国产乱子伦一区二区三区| 女性生殖器流出的白浆| 亚洲最大成人中文| 桃红色精品国产亚洲av| 亚洲全国av大片| 十分钟在线观看高清视频www| 在线观看一区二区三区| 久久午夜综合久久蜜桃| 黄色视频不卡| 欧美性长视频在线观看| 一边摸一边抽搐一进一小说| 亚洲精品国产区一区二| 国产精品秋霞免费鲁丝片| 97人妻精品一区二区三区麻豆 | 麻豆一二三区av精品| 中文字幕另类日韩欧美亚洲嫩草| 精品久久蜜臀av无| 色综合亚洲欧美另类图片| 久久人人爽av亚洲精品天堂| 91成年电影在线观看| 91九色精品人成在线观看| 色播亚洲综合网| 精品国产亚洲在线| 欧美色欧美亚洲另类二区 | 一级毛片女人18水好多| 香蕉丝袜av| 九色亚洲精品在线播放| 99香蕉大伊视频| 亚洲自拍偷在线| 国内久久婷婷六月综合欲色啪| 天堂动漫精品| 中文字幕精品免费在线观看视频| www.www免费av| 男人舔女人的私密视频| 国产亚洲欧美98| 国产又爽黄色视频| 91老司机精品| 欧美激情高清一区二区三区| 久久香蕉精品热| 热99re8久久精品国产| 中文字幕人成人乱码亚洲影| 精品久久久久久久久久免费视频| 久久中文字幕人妻熟女| 欧美成人午夜精品| 很黄的视频免费| 久久亚洲精品不卡| 亚洲精品av麻豆狂野| 精品免费久久久久久久清纯| 精品国产亚洲在线| e午夜精品久久久久久久| 精品人妻在线不人妻| 国产激情欧美一区二区| 亚洲一码二码三码区别大吗| 少妇 在线观看| 日韩大码丰满熟妇| 黑人巨大精品欧美一区二区mp4| 亚洲精品久久国产高清桃花| 国产视频一区二区在线看| a在线观看视频网站| 999久久久精品免费观看国产| 欧美成狂野欧美在线观看| 亚洲中文字幕一区二区三区有码在线看 | 日韩大码丰满熟妇| 亚洲一区中文字幕在线| 一个人观看的视频www高清免费观看 | 在线观看日韩欧美| 亚洲中文字幕一区二区三区有码在线看 | 一边摸一边抽搐一进一小说| 国产精品一区二区精品视频观看| 亚洲国产精品成人综合色| 亚洲五月色婷婷综合| 国产精品秋霞免费鲁丝片| 国产高清激情床上av| 看免费av毛片| 国产成人精品无人区| 亚洲中文字幕一区二区三区有码在线看 | 免费在线观看影片大全网站| 中出人妻视频一区二区| 老汉色av国产亚洲站长工具| 欧美日韩瑟瑟在线播放| 亚洲一卡2卡3卡4卡5卡精品中文| 中文字幕人成人乱码亚洲影| 日日夜夜操网爽| 亚洲av熟女| 亚洲一区二区三区色噜噜| www国产在线视频色| 国产精品免费一区二区三区在线| 在线观看舔阴道视频| 亚洲国产欧美一区二区综合| 大型黄色视频在线免费观看| svipshipincom国产片| 成人亚洲精品一区在线观看| 一级黄色大片毛片| 亚洲欧美激情综合另类| 黄片小视频在线播放| 国产伦一二天堂av在线观看| 久久九九热精品免费| 国产精品久久久av美女十八| 人妻丰满熟妇av一区二区三区| 亚洲精品中文字幕一二三四区| 亚洲精品中文字幕在线视频| 成熟少妇高潮喷水视频| 在线天堂中文资源库| 欧美不卡视频在线免费观看 | 禁无遮挡网站| 成人三级黄色视频| 国产99久久九九免费精品| 99国产精品一区二区蜜桃av| 18禁观看日本| 国产高清videossex| 精品久久久久久久毛片微露脸| 国产精品1区2区在线观看.| 很黄的视频免费| 久久中文看片网| 搞女人的毛片| 99精品久久久久人妻精品| 国产xxxxx性猛交| 国产精品野战在线观看| 中国美女看黄片| 午夜影院日韩av| 亚洲欧美日韩高清在线视频| 日韩av在线大香蕉| 成人国产一区最新在线观看| av天堂久久9| 欧美日韩亚洲综合一区二区三区_| 精品国产超薄肉色丝袜足j| 中出人妻视频一区二区| av有码第一页| 国产视频一区二区在线看| 久久久国产成人精品二区| 99在线人妻在线中文字幕| 欧美精品亚洲一区二区| av片东京热男人的天堂| 亚洲熟妇中文字幕五十中出| 成年版毛片免费区| 欧美色欧美亚洲另类二区 | 国产熟女xx| 制服诱惑二区| 国产主播在线观看一区二区| 亚洲国产日韩欧美精品在线观看 | 国产99久久九九免费精品| 日本 欧美在线| 91成人精品电影| 国产99久久九九免费精品| 久久天躁狠狠躁夜夜2o2o| 亚洲美女黄片视频| 啦啦啦韩国在线观看视频| 久久精品91无色码中文字幕| 精品久久久久久,| 国产精品野战在线观看| 日韩高清综合在线| 美女免费视频网站| 一区二区三区国产精品乱码| 在线观看免费日韩欧美大片| 一a级毛片在线观看| 精品国产国语对白av| 制服诱惑二区| 久久婷婷人人爽人人干人人爱 | 日韩大码丰满熟妇| 美女免费视频网站| 看黄色毛片网站| 国产精华一区二区三区| 巨乳人妻的诱惑在线观看| 国产成人精品久久二区二区免费| 成人手机av| 久热这里只有精品99| 中文字幕精品免费在线观看视频| 一进一出抽搐动态|