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

    基于馬爾科夫過(guò)程的珠江流域干旱特征研究

    2016-10-24 05:02:05肖名忠張強(qiáng)孔冬冬張正浩
    關(guān)鍵詞:珠江流域馬爾科夫站點(diǎn)

    肖名忠, 張強(qiáng), 孔冬冬, 張正浩

    (1.中山大學(xué) 水資源與環(huán)境系,廣東 廣州 510275;2.中山大學(xué) 華南地區(qū)水循環(huán)與水安全廣東省普通高校重點(diǎn)實(shí)驗(yàn)室,廣東 廣州 510275)

    ?

    基于馬爾科夫過(guò)程的珠江流域干旱特征研究

    肖名忠1,2, 張強(qiáng)1,2, 孔冬冬1,2, 張正浩1,2

    (1.中山大學(xué) 水資源與環(huán)境系,廣東 廣州 510275;2.中山大學(xué) 華南地區(qū)水循環(huán)與水安全廣東省普通高校重點(diǎn)實(shí)驗(yàn)室,廣東 廣州 510275)

    對(duì)干旱轉(zhuǎn)移概率的有效把握是干旱預(yù)報(bào)的關(guān)鍵,進(jìn)而能夠用于減輕干旱危害?;诖?,根據(jù)Copula函數(shù)導(dǎo)出了馬爾科夫鏈的轉(zhuǎn)移概率矩陣。將干旱分為了6種狀態(tài),干旱事件和干旱狀態(tài)不同,當(dāng)干旱事件的嚴(yán)重程度不小于某一干旱狀態(tài)時(shí),在本文中就被定義為具有該嚴(yán)重程度的干旱事件。為了研究珠江流域的干旱特征,對(duì)某一嚴(yán)重程度干旱事件的平均歷時(shí)、從任一干旱狀態(tài)到?jīng)]有干旱狀態(tài)的平均首達(dá)時(shí)間、以及從沒(méi)有干旱狀態(tài)到某一嚴(yán)重干旱事件的平均首達(dá)時(shí)間進(jìn)行了分析。研究結(jié)果表明:①非常極端干旱事件持續(xù)約1.5個(gè)月,輕微干旱事件持續(xù)約3個(gè)月;②從非常極端干旱狀態(tài)恢復(fù)到?jīng)]有干旱狀態(tài)需要約3.5個(gè)月,從輕微干旱狀態(tài)恢復(fù)到無(wú)旱情需要約1.7個(gè)月的時(shí)間;③一般來(lái)說(shuō),珠江流域平均0.5 a中有一次中等干旱事件,1 a中有一次嚴(yán)重干旱事件,1.5 a中有一次極端干旱事件,以及3.5 a中有一次非常極端干旱事件。此外,研究結(jié)果還表明,珠江流域東南部地區(qū)干旱風(fēng)險(xiǎn)較高,而該地區(qū)人口稠密、經(jīng)濟(jì)發(fā)達(dá),特別是其中的珠三角地區(qū),因而需要特別重視干旱風(fēng)險(xiǎn)。同時(shí),珠江流域西部地區(qū)也是干旱風(fēng)險(xiǎn)較高的區(qū)域。

    干旱;馬爾科夫鏈; Copula函數(shù);轉(zhuǎn)移概率;標(biāo)準(zhǔn)化降水蒸散發(fā)指數(shù)(SPEI)

    干旱通常定義為在某一段時(shí)間內(nèi)的持續(xù)缺水現(xiàn)象,由于持續(xù)時(shí)間長(zhǎng),干旱會(huì)導(dǎo)致水循環(huán)過(guò)程中的嚴(yán)重水量失衡,會(huì)對(duì)農(nóng)業(yè)、城市供水、公共衛(wèi)生、生態(tài)系統(tǒng)、能源以及經(jīng)濟(jì)等造成嚴(yán)重影響[1]。具體來(lái)說(shuō),干旱會(huì)對(duì)農(nóng)業(yè)、城市、供水、公共衛(wèi)生、生態(tài)系統(tǒng)、能源以及經(jīng)濟(jì)等造成嚴(yán)重影響[2-3]。據(jù)估計(jì),美國(guó)由于干旱造成的經(jīng)濟(jì)損失每年高達(dá)80億美元,遠(yuǎn)超其他氣象災(zāi)害造成的損失[4]。同時(shí),近年來(lái)隨著人口的快速增多,農(nóng)業(yè)、能源、經(jīng)濟(jì)的快速發(fā)展,對(duì)水資源的需求量正成倍增加,從而導(dǎo)致世界上很多地區(qū)幾乎每年都有缺水事件發(fā)生[1],尤其是氣候變化和水污染等因素更進(jìn)一步加劇了缺水狀況[5]。因而,干旱的成因、影響和變化規(guī)律以及風(fēng)險(xiǎn)的評(píng)估正得到各國(guó)學(xué)者的廣泛重視[6-11]。

    實(shí)際上,目前監(jiān)測(cè)干旱狀況以及評(píng)估缺水狀態(tài)的指數(shù)很多,各指數(shù)詳細(xì)的優(yōu)缺點(diǎn)可以參考Mishra和Singh寫(xiě)的綜述[1]。標(biāo)準(zhǔn)化降水指數(shù)(Standardized Precipitation Index,SPI)由于計(jì)算簡(jiǎn)單,同時(shí)能夠用不同時(shí)間尺度描述短期和長(zhǎng)期的干旱,被廣泛運(yùn)用于干旱的各個(gè)屬性的監(jiān)測(cè)[6-8, 11-13]。然而隨著全球變暖,蒸散發(fā)也增加,這將極大影響干旱的嚴(yán)重程度,因而作為SPI的一個(gè)擴(kuò)展,Vicente-Serrano等[14]提出了標(biāo)準(zhǔn)降水蒸散發(fā)指數(shù)(Standardized Precipitation Evapotranspiration Index,SPEI)。SPEI被設(shè)計(jì)成同時(shí)考慮降水和潛在蒸散發(fā)的影響監(jiān)測(cè)干旱,更能捕捉氣溫升高對(duì)干旱的影響,其在干旱的分析和評(píng)估領(lǐng)域也得到了廣泛的運(yùn)用[14-16],本文將SPEI作為干旱監(jiān)測(cè)的指標(biāo)。

    相比于洪水、地震、颶風(fēng)等極端自然災(zāi)害,干旱發(fā)生過(guò)程非常緩慢,因而在預(yù)報(bào)到位的情況下可以對(duì)干旱災(zāi)害進(jìn)行預(yù)防,從而減輕干旱的影響。同時(shí)用SPEI描述的干旱過(guò)程可近似看作為一個(gè)隨機(jī)過(guò)程,可以用隨機(jī)方法描述。Sanusi等[13]基于一階的均質(zhì)馬爾科夫過(guò)程,用1個(gè)月時(shí)間尺度的SPI對(duì)馬來(lái)西亞半島的干旱特征,如平均復(fù)發(fā)時(shí)間以及干旱的平均首達(dá)時(shí)間等進(jìn)行了分析。Paulo和Pereira[17]基于均質(zhì)和非均質(zhì)的馬爾科夫過(guò)程,對(duì)干旱狀態(tài)的轉(zhuǎn)移概率特征進(jìn)行了分析。這些研究結(jié)果發(fā)現(xiàn),隨機(jī)模型能夠很好地用于干旱演變過(guò)程的監(jiān)測(cè),從而為干旱的預(yù)防提供有用信息。

    傳統(tǒng)的基于最大似然法估計(jì)的馬爾科夫轉(zhuǎn)移概率矩陣的有效性受到了Cancelliere等[12]的質(zhì)疑。有研究發(fā)現(xiàn),由于樣本數(shù)過(guò)少,在某些情況下會(huì)導(dǎo)致計(jì)算出的轉(zhuǎn)移概率為零,這是非常不合理的。為避免傳統(tǒng)最大似然估計(jì)方法的缺陷,Cancelliere等[12]根據(jù)SPI的性質(zhì)推導(dǎo)出了計(jì)算轉(zhuǎn)移概率矩陣的理論方法。但這個(gè)方法也是有爭(zhēng)議的,Madadgar等[18]認(rèn)為該方法所基于的相互獨(dú)立以及正態(tài)分布的假設(shè)并不是經(jīng)常成立的,同時(shí)指出,根據(jù)SPI計(jì)算轉(zhuǎn)移概率矩陣的繁瑣復(fù)雜過(guò)程可以通過(guò)建立用Copula函數(shù)表示的多變量相關(guān)性模型避免。本文介紹了一種通過(guò)Copula函數(shù)計(jì)算轉(zhuǎn)移概率矩陣的理論方法。

    珠江流域地處濕潤(rùn)氣候區(qū),雨量充沛,但近年來(lái)也遭遇到了嚴(yán)重的干旱事件,例如2010年的西南5省市大旱以及2011年的南方干旱[6]。這些干旱對(duì)珠江流域的社會(huì)經(jīng)濟(jì)發(fā)展以及生態(tài)環(huán)境造成了嚴(yán)重的破壞。更重要的是,珠江流域的一些支流擔(dān)負(fù)著為珠江三角洲特大城市群以及香港供水的重任。據(jù)估計(jì),在珠江三角洲城市群生活著4 000萬(wàn)左右的居民[19],香港每年用水量的80%是由東江提供的[20]。因此,基于馬爾科夫鏈對(duì)珠江流域干旱的時(shí)空演變特征進(jìn)行分析將具有重要意義,可為決策者提供有用的信息,從而能更好地安排生產(chǎn)、生活活動(dòng),減輕干

    旱災(zāi)害目前這方面的研究在珠江流域還沒(méi)有做。

    本文的目的是利用由Copula函數(shù)計(jì)算出的馬爾科夫鏈轉(zhuǎn)移概率矩陣研究珠江流域干旱特征。研究中將干旱分為6種狀態(tài),干旱事件和干旱狀態(tài)不同,當(dāng)干旱事件的嚴(yán)重程度不小于某一干旱狀態(tài)時(shí),在本文中就被定義為具有該嚴(yán)重程度的干旱事件。對(duì)于干旱事件,了解某一嚴(yán)重程度干旱事件持續(xù)的時(shí)間、恢復(fù)到無(wú)旱情狀態(tài)的時(shí)間以及某一嚴(yán)重程度干旱事件的發(fā)生頻率將對(duì)旱災(zāi)預(yù)防及應(yīng)急管理有重大參考價(jià)值。因而,本文研究的珠江流域干旱特征是:①某一嚴(yán)重程度干旱事件的平均歷時(shí);②從任一干旱狀態(tài)到?jīng)]有干旱狀態(tài)的平均首達(dá)時(shí)間;③從沒(méi)有干旱狀態(tài)到某一嚴(yán)重干旱事件的平均首達(dá)時(shí)間。3個(gè)干旱特征的示意圖如圖1所示。

    圖1  本文分析的珠江流域干旱特征示意圖

    1 研究區(qū)域和數(shù)據(jù)

    珠江流域 (102°14′E—115°53′E;21°31′N—26°49′N) 位于中國(guó)南部,是中國(guó)河流流量第二大、流域面積第三大的流域。珠江流域面積4.573×105km2,主要由西江、北江和東江三大支流構(gòu)成.由于地處熱帶和亞熱帶季風(fēng)氣候區(qū),珠江流域多年平均氣溫在14°C到22°C,多年平均降水量1 525.1 mm,同時(shí)年內(nèi)降水分布不均,全年大約80%的降水量集中在4—9月[11]。本文的研究數(shù)據(jù)來(lái)自于珠江流域內(nèi)46個(gè)國(guó)家基本氣象站點(diǎn)從1959—2013年的日降水、最高及最低氣溫?cái)?shù)據(jù),站點(diǎn)的位置如圖2所示。

    圖2  氣象站點(diǎn)在珠江流域中的位置

    圖2中10個(gè)站點(diǎn)日降水?dāng)?shù)據(jù)有缺測(cè)值,缺測(cè)率在0.1%以下;17個(gè)站點(diǎn)日最高和最低氣溫?cái)?shù)據(jù)有缺測(cè)值,缺測(cè)率也在0.1%以下;同時(shí)對(duì)于這些有缺測(cè)數(shù)據(jù)的站點(diǎn),在圖2中也進(jìn)行了標(biāo)注。對(duì)于這些缺測(cè)的某一天的數(shù)據(jù),本文采用其他無(wú)缺測(cè)年份相同時(shí)間內(nèi)觀測(cè)值的平均值來(lái)插補(bǔ),當(dāng)實(shí)際值和均值相差過(guò)大時(shí),這種方法可能會(huì)導(dǎo)致某些有用信息的缺失,但整體影響較小。

    2 研究方法

    2.1標(biāo)準(zhǔn)化降水蒸散發(fā)指數(shù)(SPEI)

    為了加深氣溫變化對(duì)用水需求影響的了解,Vicente-Serrano等[14]提出了用SPEI指標(biāo)研究干旱。SPEI被設(shè)計(jì)成同時(shí)考慮降水(P)和潛在蒸散發(fā)(PET)的影響監(jiān)測(cè)干旱過(guò)程,其中潛在蒸散發(fā)指的是在水足夠多的情況下所能產(chǎn)生的蒸發(fā)和蒸騰總量。在SPEI計(jì)算過(guò)程中,一個(gè)改進(jìn)版的Hargreaves公式[21-23]被用來(lái)計(jì)算月潛在蒸散發(fā)。作為一個(gè)水分盈虧狀態(tài)的簡(jiǎn)單描述,月降水和潛在蒸散發(fā)的差值可以計(jì)算為Di=Pi-PETi。同時(shí)Di值在不同時(shí)間尺度下的累積值可以表示為:

    (1)

    式中:k為時(shí)間尺度;n為時(shí)間。

    和SPI值的計(jì)算類似[11, 24],SPEI的計(jì)算也是通過(guò)一個(gè)等概率轉(zhuǎn)換的方法將Dk序列轉(zhuǎn)換成一個(gè)標(biāo)準(zhǔn)正態(tài)變量。首先用不同的分布函數(shù)擬合Dk序列,選出最優(yōu)的,進(jìn)而計(jì)算非超過(guò)概率F(x),接著再將F(x)通過(guò)標(biāo)準(zhǔn)正態(tài)分布函數(shù)的反函數(shù)標(biāo)準(zhǔn)化,從而得到SPEI序列值[14]。對(duì)于一年內(nèi)不同月份的Dk序列值,分別用廣義正態(tài)分布(GNO)、廣義帕累托分布(GP)、皮爾遜三型分布(P-III)、廣義極值分布(GEV)和廣義Logistic分布(GLO)擬合,其中分布函數(shù)的參數(shù)估計(jì)采用線性矩估計(jì)[25]。根據(jù)Madadgar和Moradkhani[26]的建議,最優(yōu)分布函數(shù)是通過(guò)Kolmogorov-Smirnov(K-S)檢驗(yàn)和Akaike Information Criterion(AIC)準(zhǔn)則選取的。當(dāng)分布函數(shù)通過(guò)K-S檢驗(yàn)時(shí),AIC值最小的將選為最優(yōu)分布函數(shù)。同時(shí),和SPI的干旱分類一樣[18, 26],SPEI的干旱分類也采用美國(guó)干旱監(jiān)測(cè)中心提供的干旱分類方法,見(jiàn)表1。

    表1  基于SPEI的干旱分類[27]

    2.2Copula函數(shù)

    由于構(gòu)造聯(lián)合分布時(shí)可以不用考慮變量的邊緣分布,Copula函數(shù)已被廣泛用于多變量的相關(guān)性分析中[6, 9, 18, 26, 28]。在本文中,分別選取了來(lái)自于Archimedean(Gumbel、Clayton和Frank)、Elliptical (Gaussian和 t) 和Plackett Copula家族中的6種Copula函數(shù)來(lái)描述多變量間的相關(guān)關(guān)系[27]。對(duì)于Copula函數(shù)的擬合優(yōu)度檢驗(yàn),本文采用Cramér-von Mises 檢驗(yàn)值(Sn)選擇最優(yōu)Copula函數(shù)。Cramér-von Mises 檢驗(yàn)值定義為一個(gè)測(cè)量經(jīng)驗(yàn)Copula函數(shù)值和參數(shù)Copula函數(shù)值(CEMP和C)差距的方法[29],如下:

    (2)

    式中n為樣本量。同時(shí)該擬合優(yōu)度檢驗(yàn)方法的p值是通過(guò)自舉采樣技術(shù)(Bootstrap Sampling)得到的[26, 29-30],在本文中隨機(jī)采樣了1 000次。同樣根據(jù)Madadgar和Moradkhani[26]的建議,當(dāng)有好幾個(gè)Copula函數(shù)在95%的置信區(qū)間內(nèi)顯著時(shí),Sn值最小的Copula函數(shù)選為最優(yōu)Copula函數(shù)。

    2.3離散型馬爾科夫鏈

    離散型馬爾科夫鏈,X1、X2、…、Xn、…是一個(gè)具有馬爾科夫性質(zhì)的隨機(jī)過(guò)程,其中馬爾科夫性是指后一個(gè)狀態(tài)發(fā)生的概率只與前一個(gè)狀態(tài)有關(guān),也就是如下定義:

    (3)

    假設(shè)Pr(X1=sj|X0=si)=pij,表示的是從狀態(tài)si轉(zhuǎn)移到狀態(tài)sj的一階轉(zhuǎn)移概率,則所有狀態(tài)之間的轉(zhuǎn)移概率可以用一個(gè)轉(zhuǎn)移概率矩陣P表示如下:

    (4)

    2.4Copula函數(shù)推導(dǎo)的轉(zhuǎn)移概率

    假設(shè)X和Y分別表示的是在時(shí)刻n和n-1時(shí)的干旱事件,并分別來(lái)自于u=FX(x)和v=FY(y)的分布函數(shù)。同時(shí)干旱被分為s類,類別用Ci表示,則當(dāng)Y屬于Ci類時(shí),X屬于Cj類的條件概率可以表示如下[12]:

    (5a)

    (5b)

    另外,假設(shè)Ci和Cj類的上、下邊界分別由si1、si2、sj1、和sj2表示,且u1=FX(sj1)、u2=FX(sj2)、v1=FY(si1)、v2=FY(si2)。則當(dāng)X和Y的聯(lián)合分布用Copula函數(shù)表示的情況下,公式(5b)可以進(jìn)一步表示為:

    (6)

    需要注意的是,式(6)中的分子可以看成是對(duì)由Ci和Cj邊界構(gòu)成的方形面積的計(jì)算,同時(shí)根據(jù)Copula函數(shù)推導(dǎo)轉(zhuǎn)移概率的示意圖如圖3所示。

    圖3 基于Copula函數(shù)推導(dǎo)轉(zhuǎn)移概率的示意圖

    2.5平均首達(dá)時(shí)間

    設(shè)Tij表示狀態(tài)si首次到達(dá)狀態(tài)sj的時(shí)間,這個(gè)時(shí)間也被稱為狀態(tài)si到狀態(tài)sj的首達(dá)時(shí)間。對(duì)于Tij,它發(fā)生的概率可以定義為:

    從狀態(tài)si到狀態(tài)sj的平均首達(dá)時(shí)間Mij,也就是期望,可以計(jì)算如下:

    Mij=E[Tij=n]。

    (8)

    此外,對(duì)于干旱分析,某一程度干旱事件發(fā)生的頻率對(duì)于干旱風(fēng)險(xiǎn)評(píng)估具有重要的意義,則很有必要計(jì)算從狀態(tài)si(例如無(wú)干旱狀態(tài))到干旱事件sj+(sj+

    表示的是干旱嚴(yán)重程度不低于sj的干旱事件)的平均首達(dá)時(shí)間。和式(7)類似,狀態(tài)si到干旱事件sj+的首達(dá)時(shí)間可以表示為Tij+,從而Tij+的發(fā)生概率可以計(jì)算如下:

    式中si的嚴(yán)重程度小于sj。

    同樣,Tij+的均值(用Mij+表示)也可以類似于式(8)進(jìn)行計(jì)算。在本文中由于采用的馬爾科夫鏈?zhǔn)欠蔷|(zhì)的,Mij和Mij+不太可能通過(guò)一個(gè)解析式來(lái)進(jìn)行計(jì)算,因而本文采用隨機(jī)模擬20 000次的方法計(jì)算。

    2.6干旱平均歷時(shí)

    了解某一嚴(yán)重程度干旱事件持續(xù)的時(shí)間將會(huì)為干旱風(fēng)險(xiǎn)預(yù)防和應(yīng)急響應(yīng)提供非常有用的信息,從而減輕干旱造成的危害,因而很有必要計(jì)算某一嚴(yán)重程度干旱事件發(fā)生的平均歷時(shí)。假設(shè)干旱具有s個(gè)干旱狀態(tài),且序號(hào)小的干旱狀態(tài)表示的干旱嚴(yán)重程度低。則當(dāng)Tjj+表示干旱嚴(yán)重程度不低于第j個(gè)干旱狀態(tài)的干旱的歷時(shí)時(shí),Tjj+的發(fā)生概率可以計(jì)算如下:

    (10)

    而干旱嚴(yán)重程度不低于第j個(gè)干旱狀態(tài)的干旱的平均歷時(shí)就是Tjj+的期望,表示為E(Tjj+)。然而對(duì)于非均質(zhì)的馬爾科夫鏈,E(Tjj+)很難通過(guò)解析式計(jì)算,從而在本文中E(Tjj+)也是通過(guò)隨機(jī)模擬20 000次的方法來(lái)計(jì)算。

    3 結(jié)果和討論

    3.1分布函數(shù)和Copula函數(shù)的選擇結(jié)果

    不同時(shí)間尺度的SPEI表征的是不同類型的干旱,對(duì)應(yīng)于農(nóng)業(yè)干旱,3個(gè)月時(shí)間尺度的SPEI在本文中被用來(lái)監(jiān)測(cè)干旱狀態(tài)。本文分別采用了GNO、GP、P-III、GEV和GLO 分布函數(shù)對(duì)一年中3個(gè)月時(shí)間尺度下的各個(gè)月的Di序列進(jìn)行擬合,選擇的最優(yōu)分布函數(shù)見(jiàn)表2。從表2可以發(fā)現(xiàn),GEV、GLO和P-III分布函數(shù)對(duì)大部分站點(diǎn)各月份的Di序列擬合良好。接著根據(jù)選好的分布函數(shù),就可以基于等概率轉(zhuǎn)換方法得到SPEI值。

    表2 在對(duì)3個(gè)月尺度下每月的水分虧缺(Di)序列擬合中,具有相同最優(yōu)分布函數(shù)的站點(diǎn)個(gè)數(shù)

    另外根據(jù)3個(gè)月時(shí)間尺度SPEI值(SPEI-3)的定義,可以預(yù)期,月SPEI-3會(huì)受到前2個(gè)月的SPEI-3值的影響,同時(shí)本文也對(duì)月SPEI-3的自相關(guān)系數(shù)進(jìn)行了計(jì)算,如圖4(a)所示。從圖4(a)中可以看出,自相關(guān)系數(shù)整體上隨著通常時(shí)滯的增大而減小,而在前2個(gè)時(shí)滯的自相關(guān)性在95%的置信水平內(nèi)顯著。第二時(shí)滯的自相關(guān)性可能也受第一時(shí)滯自相關(guān)性的影響,因而本文還對(duì)月SPEI-3的部分自相關(guān)性系數(shù)進(jìn)行了計(jì)算以消除這種干擾,如圖4(b)所示。從圖4(b)中可以看出,部分自相關(guān)性在第一時(shí)滯自相關(guān)性最強(qiáng),然后隨著時(shí)滯的增加逐漸以在正負(fù)值之間震蕩的方式趨于零。因此,本文認(rèn)為月SPEI主要受前一個(gè)月值的影響,從而用一階馬爾科夫鏈表示SPEI值的狀態(tài)轉(zhuǎn)移概率,同時(shí)連續(xù)2個(gè)月的SPEI-3值間的相關(guān)性用Copula函數(shù)表示。根據(jù)2.2部分介紹的Copula擬合優(yōu)度檢驗(yàn)方法,本文選擇的最優(yōu)Copula函數(shù)見(jiàn)表3,在用Copula函數(shù)對(duì)連續(xù)2個(gè)月相關(guān)關(guān)系的描述中,具有相同Copula函數(shù)的站點(diǎn)個(gè)數(shù),其中連續(xù)2個(gè)月中的第一個(gè)月份顯示在表3中。從表3中可以看出Gumbel、Gaussian和t Copula函數(shù)對(duì)連續(xù)2個(gè)月SPEI-3相關(guān)性的描述表現(xiàn)不錯(cuò)。此外,由于月SPEI-3序列是標(biāo)準(zhǔn)正態(tài)變量,因而本文中Copula函數(shù)的邊緣分布為正態(tài)分布。

    圖4 各站點(diǎn)3個(gè)月時(shí)間尺度SPEI的箱型圖(紅點(diǎn)線表示的是95%置信區(qū)間的閾值)

    Copula函數(shù)1月2月3月4月5月6月7月8月9月10月11月12月Gumbel1620117313191278159Clayton3—79131—2—811Frank44261332154103Gaussian9714141614121713121120t10121071313912613610Plackett3213—2312132

    3.2干旱特征的空間分布

    通常馬爾科夫鏈的轉(zhuǎn)移概率矩陣是根據(jù)觀測(cè)到的轉(zhuǎn)移樣本以最大似然法進(jìn)行估算的[13, 17],同樣,本文也根據(jù)最大似然法對(duì)各站點(diǎn)的轉(zhuǎn)移概率矩陣進(jìn)行了估算。然而由于樣本量有限,在某些情況下會(huì)得到轉(zhuǎn)移概率為零的情況,見(jiàn)表4。在表1中,本文定義了6種干旱狀態(tài)(無(wú)干旱、D0、D1、D2、D3和D4),依次表示干旱嚴(yán)重程度從無(wú)干旱到非常極端的狀態(tài)。從表4中可以看出很多站點(diǎn)從D0到D4、D4到D0和D4到無(wú)干旱狀態(tài)的轉(zhuǎn)移概率都為零,這明顯是不符合實(shí)際的。因而為了避免這種缺陷,本文采用在前文中選好的最優(yōu)Copula函數(shù)從理論上計(jì)算各站點(diǎn)各月份的轉(zhuǎn)移概率矩陣。結(jié)果發(fā)現(xiàn),對(duì)于所有站點(diǎn),每個(gè)月的狀態(tài)轉(zhuǎn)移概率都是變化的(圖5),因而本文采用的馬爾科夫鏈?zhǔn)欠蔷|(zhì)的。接著基于這非均質(zhì)馬爾科夫鏈,對(duì)某一嚴(yán)重程度干旱事件的平均歷時(shí)、從任一干旱狀態(tài)到?jīng)]有干旱狀態(tài)的平均首達(dá)時(shí)間以及從沒(méi)有干旱狀態(tài)到某一嚴(yán)重干旱事件的平均首達(dá)時(shí)間分別進(jìn)行了計(jì)算。

    表4最大似然法估計(jì)的轉(zhuǎn)移概率中概率值為0的站點(diǎn)個(gè)數(shù)

    D0D1D2D3D4無(wú)干旱D03065260D1000060D2301080D3300010D425431417無(wú)干旱000170

    圖5 珠江流域各站點(diǎn)從任何干旱狀態(tài)到D4狀態(tài)和無(wú)旱情狀態(tài)的各月轉(zhuǎn)移概率變化情況

    通過(guò)反距離權(quán)重插值方法,分析了干旱平均歷時(shí)的空間分布情況,如圖6所示。從圖6中可以看出,珠江流域非常極端干旱的平均歷時(shí)約1.5個(gè)月,輕微干旱的平均歷時(shí)約3個(gè)月,這些結(jié)果將會(huì)為珠江流域的干旱預(yù)防和應(yīng)急管理提供有用的信息。同時(shí),從圖6中可以看到珠江流域東南部干旱的平均歷時(shí)更長(zhǎng),尤其是極端干旱。圖7給出了從任一干旱狀態(tài)到無(wú)干旱狀態(tài)的平均首達(dá)時(shí)間在珠江流域的空間分布,從圖中也發(fā)現(xiàn),珠江流域東南部從干旱狀態(tài)到無(wú)旱情狀態(tài)的平均首達(dá)時(shí)間更長(zhǎng)。這些說(shuō)明珠江流域東南部遭受干旱的概率更高,這也與我們之前的文章得到的珠江流域東南部干旱風(fēng)險(xiǎn)更高的結(jié)果一致[6, 11]。珠江流域東南部作為一個(gè)人口密集、經(jīng)濟(jì)發(fā)達(dá)的區(qū)域,高發(fā)的干旱風(fēng)險(xiǎn)將會(huì)對(duì)該區(qū)域社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展造成嚴(yán)重的影響,尤其是對(duì)珠三角地區(qū),因而需引起高度重視。類似的,從圖6—7中可以看出珠江流域的西部地區(qū)也是一個(gè)干旱風(fēng)險(xiǎn)高的地區(qū)。

    圖6 某一嚴(yán)重程度干旱的平均歷時(shí)在珠江流域的空間分布圖(圖中顏色條的單位是月)

    圖7 從任一干旱狀態(tài)到達(dá)無(wú)干旱狀態(tài)的平均首達(dá)時(shí)間在珠江流域的空間分布圖(圖中顏色條的單位是月)

    從無(wú)干旱狀態(tài)到某一嚴(yán)重干旱事件的平均首達(dá)時(shí)間的空間分布情況如圖8所示。從圖8中可看出,整體上看,從無(wú)干旱狀態(tài)到某一嚴(yán)重干旱事件在珠江流域的平均首達(dá)時(shí)間對(duì)于中等干旱是0.5 a、嚴(yán)重干旱是1 a、極端干旱是1.5 a、非常極端干旱是3.5 a。同時(shí),珠江流域上游地區(qū)干旱的發(fā)生率更高,這進(jìn)一步增加了珠江流域西部地區(qū)的干旱風(fēng)險(xiǎn)。另外,珠江流域東南部和西部的干旱一旦發(fā)生,它們將持續(xù)更長(zhǎng)時(shí)間(圖6—8)。綜合上述所有干旱特征信息,可以更好地了解珠江流域各地區(qū)的干旱變化情況,從而為水資源的管理和決策提供參考信息。

    圖8 從無(wú)干旱狀態(tài)到達(dá)某一嚴(yán)重程度干旱的平均首達(dá)時(shí)間在珠江流域的空間分布圖(圖中顏色條的單位是月)

    4 結(jié) 語(yǔ)

    很多計(jì)算發(fā)現(xiàn)由于樣本量有限,在某些情況下會(huì)導(dǎo)致根據(jù)最大似然法計(jì)算出的轉(zhuǎn)移概率為零,這是非常不合理的。為了避免傳統(tǒng)最大似然估計(jì)方法的缺陷,本文介紹了一種用Copula函數(shù)計(jì)算各月轉(zhuǎn)移概率矩陣的方法,同時(shí)由于計(jì)算出的轉(zhuǎn)移概率是隨時(shí)間變化的,因而采用的馬爾科夫鏈?zhǔn)欠蔷|(zhì)的。接著,根據(jù)計(jì)算出的非均質(zhì)馬爾科夫鏈,以珠江流域?yàn)檠芯繉?duì)象,對(duì)某一嚴(yán)重程度干旱事件的平均歷時(shí)、從任一干旱狀態(tài)到?jīng)]有干旱狀態(tài)的平均首達(dá)時(shí)間以及從沒(méi)有干旱狀態(tài)到某一嚴(yán)重干旱事件的平均首達(dá)時(shí)間3個(gè)干旱特征分別進(jìn)行了分析。結(jié)果發(fā)現(xiàn):①非常極端干旱事件持續(xù)約1.5個(gè)月,輕微干旱事件持續(xù)約3個(gè)月;②從非常極端干旱狀態(tài)恢復(fù)到無(wú)干旱狀態(tài)需要約3.5個(gè)月時(shí)間,而從輕微干旱狀態(tài)恢復(fù)到無(wú)旱情需要約1.7個(gè)月的時(shí)間;③一般來(lái)說(shuō),珠江流域平均0.5 a中發(fā)生一次中等干旱事件,1 a中發(fā)生一次嚴(yán)重干旱事件,1.5 a中發(fā)生一次極端干旱事件,以及3.5 a中發(fā)生一次非常極端干旱事件。這些信息將會(huì)對(duì)珠江流域的干旱災(zāi)害預(yù)防及應(yīng)急管理提供非常有用的信息。此外,研究結(jié)果還表明,珠江流域東南部地區(qū)干旱風(fēng)險(xiǎn)較高,而該地區(qū)人口稠密、經(jīng)濟(jì)發(fā)達(dá),特別是其中的珠三角地區(qū),因而需要加強(qiáng)對(duì)該地區(qū)的干旱風(fēng)險(xiǎn)評(píng)估和管理。同時(shí),珠江流域西部地區(qū)也是干旱風(fēng)險(xiǎn)較高的區(qū)域。

    [1]Mishra A K,Singh V P.A review of drought concepts[J].Journal of Hydrology,2010,391(2):202-216.

    [2]Devereux S.The impact of droughts and floods on food security and policy options to alleviate negative effects[J].Agricultural Economics,2007,37(1):47-58.

    [3]Zhang Q,Sun P,Li J F,et al.Spatiotemporal properties of droughts and related impacts on agriculture in Xinjiang, China[J].International Journal of Climatology,2015,35(7):1254-1266.

    [4]Wilhite D A.Drought as a natural hazard:concepts and definitions[M]//Drought:A Global Assessment[M].London:Routledge Publishers,2000:3-18.

    [5]V?r?smarty C J,Green P,Salisbury J,et al.Global water resources:vulnerability from climate change and population growth[J].Science,2000,289:284-288.

    [6]Zhang Q,Xiao M Z,Singh V P,et al.Regionalization and spatial changing properties of droughts across the Pearl River basin,China[J].Journal of Hydrology,2012,472-473:355-366.

    [7]Kao S C,Govindaraju R S.A copula-based joint deficit index for droughts[J].Journal of Hydrology,2010,380(1-2):121-134.

    [8]Shiau J.Fitting drought duration and severity with two-dimensional copulas[J].Water Resources Management,2006,20(5):795-815.

    [9]Zhang Q,Xiao M Z,Singh V,et al.Copula-based risk evaluation of hydrological droughts in the East River basin,China[J].Stochastic Environmental Research and Risk Assessment,2013,27(6):1397-1406.

    [10]Loukas A,Vasiliades L.Probabilistic analysis of drought spatiotemporal characteristics in Thessaly region,Greece[J].Natural Hazards and Earth System Science,2004,4:719-731.

    [11]肖名忠,張強(qiáng),陳曉宏.基于多變量概率分析的珠江流域干旱特征研究 [J].地理學(xué)報(bào),2012,67(1):83-92.

    [12]Cancelliere A,Mauro G D,Bonaccorso B,et al.Drought forecasting using the standardized precipitation index [J].Water Resources Management,2007,21(5):801-819.

    [13]Sanusi W,Jemain A,Zin W,et al.The drought characteristics using the first-order homogeneous markov chain of monthly rainfall data in peninsular malaysia[J].Water Resources Management,2015,29(5):1523-1539.

    [14]Vicente-Serrano S M,Beguer A S,López-moreno J I.A multiscalar drought index sensitive to global warming:The standardized precipitation evapotranspiration index[J].Journal of Climate,2010,23(7):1696-1718.

    [15]Yu M,Li Q,Hayes M J,et al.Are droughts becoming more frequent or severe in China based on the standardized precipitation evapotranspiration index:1951—2010? [J].International Journal of Climatology,2014,34(3):545-558.

    [16]Zhang Q,Qi T Y,Singh V P,et al.Regional frequency analysis of droughts in China:multivariate perspective[J].Water Resources Management,2015,29(6):1767-1787.

    [17]Paulo A,Pereira L S.Prediction of SPI drought class transitions using markov chains[J].Water Resources Management,2007,21(10):1813-1827.

    [18]Madadgar S,Moradkhani H.A bayesian framework for probabilistic seasonal drought forecasting[J].Journal of Hydrometeorology,2013,14(6):1685-1705.

    [19]Yang T,Shao Q,Hao Z C,et al.Regional frequency analysis and spatio-temporal pattern characterization of rainfall extremes in the Pearl River Basin,China[J].Journal of Hydrology,2010,380(3-4):386-405.

    [20]Wong J S,Zhang Q,Chen Y D.Statistical modeling of daily urban water consumption in Hong Kong:trend,changing patterns,and forecast[J].Water Resources Research,2010,46(3):W03506.

    [21]Droogers P,Allen R.Estimating reference evapotranspiration under inaccurate data conditions[J].Irrigation and Drainage Systems,2002,16(1):33-45.

    [22]Hargreaves G H.Defining and using reference evapotranspiration[J].Journal of Irrigation and Drainage Engineering,1994,120(6):1132-1139.

    [23]Beguer A S,Vicente-Serrano S M.SPEI:Calculation of the Standardised Precipitation-Evapotranspiration Index[R/OL].[2015-01-20].http://CRAN.R-project.org/package=SPEI.

    [24]Mckee T B,Doesken N J,KLEIST J.The relationship of drought frequency and duration to time scales[C]//Proceedings of the Proc Eighth Conf on Applied Climatology.Amer.Meteor.Soc.Anaheim,CA,F:1993.

    [25]Hosking J R M.L-moments:analysis and estimation of distributions using linear combinations of order statistics[J].Journal of the Royal Statistical Society Series B (Methodological),1990,52(1):105-124.

    [26]Madadgar S,Moradkhani H.Spatio-temporal drought forecasting within Bayesian networks[J].Journal of Hydrology,2014,512:134-146.

    [27]Chen Y D,Zhang Q,Xiao M Z,et al.Probabilistic forecasting of seasonal droughts in the Pearl River basin,China[J/OL].Stochastic Environmental Research and Risk Assessment,2015:1-10(2015-11-20)[2015-10-24].http:∥link.springer.com/journal/477.DOI:10.1007/s00477-015-1174-6.

    [28]Nelsen R B.An Introduction to Copulas[M].New York:Springer Verlag,2006.

    [29]Genest C,R Millard B,Beaudoin D.Goodness-of-fit tests for copulas:a review and a power study[J].Insurance:Mathematics and Economics,2009,44(2):199-213.

    [30]Hofert M,Kojadinovic I,Maechler M,et al.Copula:Multivariate Dependence with Copulas[R/OL].[2015-10-24].http:∥CRAN.R-project.org/package=copula.

    (責(zé)任編輯:蔡洪濤)

    Research on Drought Behaviors in the Pearl River Basin Based on Markov Chain

    XIAO Mingzhong1, 2, ZHANG Qiang1, 2, KONG Dongdong1, 2, ZHANG Zhenghao1, 2

    (1. Department of Water Resources and Environment, Sun Yat-sen University, Guangzhou 510275, China;2. Key Laboratory of Water Cycle and Water Security in Southern China of Guangdong High Education Institute,Sun Yat-sen University, Guangzhou 510275, China)

    Deep understanding of transition probability behaviors of droughts is the key to forecast of droughts, and mitigation of drought hazards. In this case, the transition probability matrix of the Markov chain was derived by the Copula functions. Six drought status were classified in the study, drought event was different from drought status, and a drought event with the severity not less than a certain drought status had been considered as a certain severe drought. In order to investigate the drought behaviors in the Pearl River basin, the mean duration of a certain severe drought, the mean first passage time from any drought status to no drought status, and the mean first passage time from no drought status to a certain severe drought were analyzed in this study. The results of this study indicated that: (1) the exceptional drought event lasted about 1.5 months, and the mild drought event lasted about 3 months; (2) it took about 3.5 months from the exceptional drought status to no drought status, while it took about 1.7 months from the mild drought status to no drought status; (3) generally, and there was a moderate drought event average in 0.5 years, a severe drought event in 1 year, an extreme drought event in 1.5 years, and an exceptional drought event in 3.5 years. Further, the results of this study indicated higher risk of drought in the southeast part of the Pearl River basin, and further works needed to be done to prevent or mitigate the drought hazards, as these regions were densely populated and economically developed, especially in the Pearl River delta. In addition, the west part of the Pearl River basin had also been identified as a region with a higher risk of drought.

    drought; Markov chain; copula functions; transition probability; standardized precipitation evapotranspiration index (SPEI)

    2015-12-20

    國(guó)家杰出青年科學(xué)基金項(xiàng)目(51425903);香港特別行政區(qū)研究資助局項(xiàng)目(CUHK441313);中央高?;究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金聯(lián)合資助。

    肖名忠(1990—),男,江西吉安人,博士,主要從事氣象水文極值分析方面的研究。E-mail:xmingzh@mail2.sysu.edu.cn。

    張強(qiáng)(1974—),男,山東沂水人,教授,博士生導(dǎo)師,博士,主要從事流域氣象水文學(xué)、旱澇災(zāi)害機(jī)理、流域地表水文過(guò)程及其對(duì)氣候變化的響應(yīng)機(jī)制與機(jī)理以及流域生態(tài)需水等領(lǐng)域的研究。E-mail:zhangq68@mail.sysu.edu.cn。

    10.3969/j.issn.1002-5634.2016.01.003

    TV882.4;P333

    A

    1002-5634(2016)01-0015-09

    猜你喜歡
    珠江流域馬爾科夫站點(diǎn)
    珠江流域自然水面綠色生態(tài)增殖技術(shù)探討
    基于CMIP6的珠江流域未來(lái)干旱時(shí)空變化
    基于疊加馬爾科夫鏈的邊坡位移預(yù)測(cè)研究
    基于改進(jìn)的灰色-馬爾科夫模型在風(fēng)機(jī)沉降中的應(yīng)用
    基于Web站點(diǎn)的SQL注入分析與防范
    電子制作(2019年14期)2019-08-20 05:43:42
    2017~2018年冬季西北地區(qū)某站點(diǎn)流感流行特征分析
    我國(guó)珠江流域水文化遺產(chǎn)數(shù)字化保護(hù)策略
    首屆歐洲自行車共享站點(diǎn)協(xié)商會(huì)召開(kāi)
    怕被人認(rèn)出
    馬爾科夫鏈在教學(xué)評(píng)價(jià)中的應(yīng)用
    91成人精品电影| 欧美日韩视频精品一区| 丝袜人妻中文字幕| 欧美日韩福利视频一区二区| 久久久久国产一级毛片高清牌| 两性夫妻黄色片| 午夜91福利影院| 99久久综合免费| 国产99久久九九免费精品| 免费av中文字幕在线| 午夜免费鲁丝| 一本久久精品| 国产精品久久久久成人av| 欧美精品亚洲一区二区| 天堂俺去俺来也www色官网| 国产精品一二三区在线看| 国产日韩欧美在线精品| 久久天堂一区二区三区四区| 大香蕉久久网| 人人妻人人爽人人添夜夜欢视频| 最近中文字幕2019免费版| 亚洲精品美女久久av网站| 久久久国产精品麻豆| 久久精品人人爽人人爽视色| 亚洲人成77777在线视频| 看十八女毛片水多多多| 中文字幕另类日韩欧美亚洲嫩草| 高清黄色对白视频在线免费看| 一级毛片黄色毛片免费观看视频| 午夜福利免费观看在线| 亚洲人成77777在线视频| 成人亚洲欧美一区二区av| 亚洲成人av在线免费| 免费观看人在逋| 久久久久久久久久久久大奶| 日本一区二区免费在线视频| av网站在线播放免费| 一本大道久久a久久精品| 十八禁高潮呻吟视频| 看免费成人av毛片| 成人黄色视频免费在线看| 久久ye,这里只有精品| 曰老女人黄片| 另类亚洲欧美激情| 在线看a的网站| 日韩中文字幕视频在线看片| 美女国产高潮福利片在线看| 精品一品国产午夜福利视频| 亚洲精品日韩在线中文字幕| 久久久久久久久免费视频了| 在线观看www视频免费| 亚洲人成网站在线观看播放| 欧美亚洲 丝袜 人妻 在线| 国产视频首页在线观看| 少妇人妻精品综合一区二区| 免费在线观看完整版高清| 视频在线观看一区二区三区| 黄片播放在线免费| 亚洲精品国产区一区二| 日本黄色日本黄色录像| 看非洲黑人一级黄片| 黄色毛片三级朝国网站| 精品亚洲乱码少妇综合久久| 大陆偷拍与自拍| 国内精品久久久久久久电影| 欧美国产精品va在线观看不卡| 波多野结衣av一区二区av| 美女高潮喷水抽搐中文字幕| 国产又爽黄色视频| 亚洲一区中文字幕在线| 妹子高潮喷水视频| 成人av一区二区三区在线看| 一级a爱片免费观看的视频| 嫩草影视91久久| 精品久久久久久,| 午夜福利一区二区在线看| 免费高清视频大片| 成在线人永久免费视频| 一边摸一边抽搐一进一小说| 一区在线观看完整版| 日本撒尿小便嘘嘘汇集6| 亚洲精品国产一区二区精华液| 淫妇啪啪啪对白视频| av电影中文网址| 久久香蕉精品热| 午夜影院日韩av| www.999成人在线观看| 最近最新中文字幕大全免费视频| 中出人妻视频一区二区| 国产精品98久久久久久宅男小说| 亚洲久久久国产精品| 久9热在线精品视频| 中文字幕av电影在线播放| 大型黄色视频在线免费观看| 精品国产乱码久久久久久男人| 最新在线观看一区二区三区| 女性被躁到高潮视频| 一a级毛片在线观看| a级毛片在线看网站| 最新美女视频免费是黄的| 久久精品91蜜桃| 国产欧美日韩精品亚洲av| 69精品国产乱码久久久| 精品卡一卡二卡四卡免费| av中文乱码字幕在线| 91字幕亚洲| av有码第一页| 国产单亲对白刺激| 国产乱人伦免费视频| 两个人免费观看高清视频| 最近最新免费中文字幕在线| 精品久久久久久久久久免费视频| 国产精品久久久久久亚洲av鲁大| 久久久久久久久久久久大奶| 黄片播放在线免费| 人人妻,人人澡人人爽秒播| 俄罗斯特黄特色一大片| 人妻久久中文字幕网| 日韩精品中文字幕看吧| 亚洲av美国av| 欧美人与性动交α欧美精品济南到| 免费一级毛片在线播放高清视频 | 88av欧美| 欧美黄色淫秽网站| 一边摸一边抽搐一进一出视频| 黄色女人牲交| 日本一区二区免费在线视频| 天天一区二区日本电影三级 | 我的亚洲天堂| 无遮挡黄片免费观看| 满18在线观看网站| 两个人看的免费小视频| 嫩草影视91久久| 日韩欧美免费精品| 国内久久婷婷六月综合欲色啪| 亚洲人成网站在线播放欧美日韩| АⅤ资源中文在线天堂| 少妇的丰满在线观看| av免费在线观看网站| 男女床上黄色一级片免费看| 久久久国产成人精品二区| 亚洲av五月六月丁香网| 国产欧美日韩精品亚洲av| 人人妻人人澡人人看| 亚洲专区字幕在线| 18禁国产床啪视频网站| 国产精品秋霞免费鲁丝片| 国产免费男女视频| 黄色女人牲交| 99久久国产精品久久久| 激情在线观看视频在线高清| 国产一区二区在线av高清观看| 性少妇av在线| 久久香蕉精品热| 天堂√8在线中文| 18禁黄网站禁片午夜丰满| 又黄又粗又硬又大视频| 91麻豆精品激情在线观看国产| 久久香蕉激情| 老司机在亚洲福利影院| 欧美日本中文国产一区发布| 国产区一区二久久| 国产视频一区二区在线看| 久久精品影院6| 日韩一卡2卡3卡4卡2021年| 中文字幕精品免费在线观看视频| 亚洲av成人av| 成熟少妇高潮喷水视频| 老熟妇仑乱视频hdxx| 久久人妻熟女aⅴ| 国产人伦9x9x在线观看| 999久久久国产精品视频| 丝袜美腿诱惑在线| 亚洲欧洲精品一区二区精品久久久| 99久久久亚洲精品蜜臀av| 叶爱在线成人免费视频播放| 久久久久精品国产欧美久久久| 999久久久国产精品视频| 亚洲精品中文字幕一二三四区| 日韩欧美国产在线观看| 久久久久亚洲av毛片大全| cao死你这个sao货| 久久国产亚洲av麻豆专区| 亚洲国产精品sss在线观看| 1024视频免费在线观看| 看黄色毛片网站| 国产又色又爽无遮挡免费看| 久久午夜亚洲精品久久| 亚洲天堂国产精品一区在线| 老司机靠b影院| 精品国产乱码久久久久久男人| 中文字幕最新亚洲高清| 免费看a级黄色片| 亚洲精品久久国产高清桃花| 亚洲午夜理论影院| 欧美午夜高清在线| 日本一区二区免费在线视频| 久热爱精品视频在线9| av有码第一页| 校园春色视频在线观看| 在线观看免费视频日本深夜| 一区二区日韩欧美中文字幕| 久久久精品欧美日韩精品| 亚洲人成电影观看| 成人亚洲精品av一区二区| 亚洲国产欧美网| 久久久久久人人人人人| 1024香蕉在线观看| 久久精品成人免费网站| 欧美中文综合在线视频| 欧美国产精品va在线观看不卡| 女人高潮潮喷娇喘18禁视频| 天堂动漫精品| 成年人黄色毛片网站| 久久国产精品人妻蜜桃| 黑人巨大精品欧美一区二区mp4| 精品一品国产午夜福利视频| 免费人成视频x8x8入口观看| 大码成人一级视频| 亚洲自偷自拍图片 自拍| 国产xxxxx性猛交| av视频在线观看入口| 亚洲欧美日韩高清在线视频| 国产不卡一卡二| 欧美日韩一级在线毛片| 侵犯人妻中文字幕一二三四区| 久久久久久久午夜电影| 国产aⅴ精品一区二区三区波| 日本五十路高清| 嫁个100分男人电影在线观看| 国产亚洲精品一区二区www| 久久中文看片网| 久久精品国产99精品国产亚洲性色 | 欧美绝顶高潮抽搐喷水| 色综合站精品国产| 精品第一国产精品| 国产又爽黄色视频| 免费在线观看亚洲国产| 黄色 视频免费看| 大型av网站在线播放| 可以在线观看毛片的网站| 国产高清激情床上av| 老司机福利观看| 精品电影一区二区在线| 在线观看午夜福利视频| 色婷婷久久久亚洲欧美| 亚洲av电影在线进入| 日本三级黄在线观看| 丝袜人妻中文字幕| 精品国产乱子伦一区二区三区| 99re在线观看精品视频| 欧美色欧美亚洲另类二区 | 在线观看免费午夜福利视频| 黄网站色视频无遮挡免费观看| 久久中文字幕人妻熟女| 91精品国产国语对白视频| 欧美午夜高清在线| 欧美另类亚洲清纯唯美| 欧美av亚洲av综合av国产av| 欧美人与性动交α欧美精品济南到| 午夜激情av网站| 国产精品自产拍在线观看55亚洲| 麻豆一二三区av精品| 少妇被粗大的猛进出69影院| 色综合欧美亚洲国产小说| 日本免费一区二区三区高清不卡 | 午夜福利高清视频| 这个男人来自地球电影免费观看| 9热在线视频观看99| 久久久精品欧美日韩精品| 淫妇啪啪啪对白视频| 精品一品国产午夜福利视频| 色综合欧美亚洲国产小说| 波多野结衣高清无吗| 青草久久国产| 精品少妇一区二区三区视频日本电影| 男女下面插进去视频免费观看| 久久久久久久久中文| 国产精品免费一区二区三区在线| 国产精品永久免费网站| avwww免费| 日韩一卡2卡3卡4卡2021年| 亚洲一区二区三区色噜噜| 久久久久久久久久久久大奶| 又黄又爽又免费观看的视频| а√天堂www在线а√下载| 麻豆一二三区av精品| 岛国视频午夜一区免费看| 国产成人啪精品午夜网站| 涩涩av久久男人的天堂| 亚洲人成电影免费在线| 黄色 视频免费看| 18美女黄网站色大片免费观看| 老司机深夜福利视频在线观看| 亚洲av日韩精品久久久久久密| 91老司机精品| 国产精品,欧美在线| 国产亚洲欧美精品永久| 欧美av亚洲av综合av国产av| 国产精品亚洲av一区麻豆| 精品电影一区二区在线| 亚洲人成电影免费在线| 每晚都被弄得嗷嗷叫到高潮| 午夜精品在线福利| 亚洲精品在线美女| 嫁个100分男人电影在线观看| 亚洲色图av天堂| 免费观看人在逋| 在线国产一区二区在线| 国产三级黄色录像| 国产成人免费无遮挡视频| 两个人免费观看高清视频| 亚洲第一欧美日韩一区二区三区| 国产精品久久久av美女十八| 精品一区二区三区四区五区乱码| 日韩免费av在线播放| 在线观看免费日韩欧美大片| 激情视频va一区二区三区| 亚洲av第一区精品v没综合| 亚洲午夜理论影院| 天天添夜夜摸| 黑人操中国人逼视频| 国产黄a三级三级三级人| 欧美大码av| 久热这里只有精品99| 欧美成人免费av一区二区三区| 99国产极品粉嫩在线观看| av福利片在线| 亚洲专区中文字幕在线| 热99re8久久精品国产| 少妇粗大呻吟视频| 69精品国产乱码久久久| 国产精品综合久久久久久久免费 | 免费无遮挡裸体视频| 黄频高清免费视频| 欧美国产精品va在线观看不卡| 亚洲九九香蕉| 亚洲一区二区三区不卡视频| 久久人妻福利社区极品人妻图片| 又黄又粗又硬又大视频| 亚洲欧美一区二区三区黑人| www日本在线高清视频| 99国产综合亚洲精品| 麻豆一二三区av精品| 19禁男女啪啪无遮挡网站| 九色国产91popny在线| 国产av一区二区精品久久| 欧美成人免费av一区二区三区| 亚洲欧美日韩无卡精品| 欧美成人免费av一区二区三区| 亚洲国产毛片av蜜桃av| 狂野欧美激情性xxxx| 色在线成人网| 亚洲成国产人片在线观看| 嫩草影视91久久| 18禁黄网站禁片午夜丰满| 老熟妇仑乱视频hdxx| 午夜精品久久久久久毛片777| 不卡一级毛片| 精品熟女少妇八av免费久了| 韩国精品一区二区三区| e午夜精品久久久久久久| ponron亚洲| 亚洲aⅴ乱码一区二区在线播放 | 91成年电影在线观看| 亚洲第一欧美日韩一区二区三区| 夜夜看夜夜爽夜夜摸| 1024香蕉在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 国产亚洲精品综合一区在线观看 | 亚洲免费av在线视频| 国产精品一区二区精品视频观看| 国产黄a三级三级三级人| 大码成人一级视频| 动漫黄色视频在线观看| 天天添夜夜摸| 欧美大码av| 午夜免费成人在线视频| 1024视频免费在线观看| 亚洲av片天天在线观看| 久久婷婷人人爽人人干人人爱 | 九色亚洲精品在线播放| 制服丝袜大香蕉在线| 咕卡用的链子| 国产精品自产拍在线观看55亚洲| 香蕉久久夜色| 丝袜美腿诱惑在线| 免费在线观看视频国产中文字幕亚洲| 制服丝袜大香蕉在线| 国产1区2区3区精品| 国产三级黄色录像| 亚洲av片天天在线观看| 久久中文字幕一级| 午夜福利影视在线免费观看| 久久人妻熟女aⅴ| 亚洲专区字幕在线| 男女下面插进去视频免费观看| 不卡av一区二区三区| 免费女性裸体啪啪无遮挡网站| 国产成+人综合+亚洲专区| 成在线人永久免费视频| 欧美绝顶高潮抽搐喷水| 亚洲aⅴ乱码一区二区在线播放 | 国产精品一区二区在线不卡| 午夜福利,免费看| 人人澡人人妻人| 日本精品一区二区三区蜜桃| 亚洲精品av麻豆狂野| 黄色视频不卡| 免费一级毛片在线播放高清视频 | 国产伦一二天堂av在线观看| 18禁美女被吸乳视频| 欧美性长视频在线观看| 亚洲精华国产精华精| 一二三四在线观看免费中文在| 午夜福利在线观看吧| 亚洲精品一卡2卡三卡4卡5卡| 久久天躁狠狠躁夜夜2o2o| 午夜激情av网站| 18美女黄网站色大片免费观看| 在线观看免费视频日本深夜| 一区二区三区精品91| 99riav亚洲国产免费| 国产亚洲精品久久久久久毛片| 高清毛片免费观看视频网站| 国产成人欧美| 国产成人av教育| 乱人伦中国视频| 久久久久久亚洲精品国产蜜桃av| 国产在线观看jvid| 精品久久久久久成人av| 淫秽高清视频在线观看| 亚洲成av片中文字幕在线观看| 欧美日韩亚洲国产一区二区在线观看| 一a级毛片在线观看| 欧美色欧美亚洲另类二区 | 欧美中文日本在线观看视频| 日本 av在线| 天堂影院成人在线观看| 亚洲精品国产精品久久久不卡| 欧美日本视频| 夜夜爽天天搞| 麻豆久久精品国产亚洲av| 三级毛片av免费| 亚洲美女黄片视频| 日本在线视频免费播放| 天堂√8在线中文| 99久久99久久久精品蜜桃| 黄频高清免费视频| 亚洲自偷自拍图片 自拍| 一级片免费观看大全| 亚洲熟女毛片儿| 亚洲一卡2卡3卡4卡5卡精品中文| 村上凉子中文字幕在线| 欧洲精品卡2卡3卡4卡5卡区| 免费高清视频大片| 久久精品国产综合久久久| 狠狠狠狠99中文字幕| 久热爱精品视频在线9| 好男人在线观看高清免费视频 | 国产精品av久久久久免费| 在线十欧美十亚洲十日本专区| 日本 欧美在线| 久久精品91蜜桃| 人人妻人人澡人人看| 亚洲欧美日韩无卡精品| 天堂影院成人在线观看| www.999成人在线观看| 50天的宝宝边吃奶边哭怎么回事| 黄色女人牲交| 啦啦啦 在线观看视频| 母亲3免费完整高清在线观看| 在线十欧美十亚洲十日本专区| 亚洲九九香蕉| 波多野结衣巨乳人妻| 国产三级黄色录像| 午夜亚洲福利在线播放| 后天国语完整版免费观看| 亚洲熟女毛片儿| www.www免费av| 国产精品免费视频内射| 69av精品久久久久久| 欧美性长视频在线观看| 亚洲黑人精品在线| 亚洲最大成人中文| 国产精品自产拍在线观看55亚洲| 波多野结衣巨乳人妻| 99在线人妻在线中文字幕| 人人澡人人妻人| 99热只有精品国产| 老司机福利观看| 最好的美女福利视频网| 女警被强在线播放| 一区二区三区高清视频在线| 欧美成狂野欧美在线观看| 黄网站色视频无遮挡免费观看| 1024视频免费在线观看| 亚洲欧美精品综合一区二区三区| 久久久久久久精品吃奶| 日韩国内少妇激情av| 丝袜美腿诱惑在线| 欧美日本视频| 丝袜美腿诱惑在线| 亚洲一区二区三区不卡视频| 极品教师在线免费播放| 精品国产国语对白av| or卡值多少钱| 在线观看舔阴道视频| 国产av精品麻豆| 久久精品亚洲精品国产色婷小说| 黄色视频,在线免费观看| 亚洲一区二区三区不卡视频| 中文字幕另类日韩欧美亚洲嫩草| 一级毛片精品| 国产成人精品在线电影| 黑人操中国人逼视频| 中文字幕另类日韩欧美亚洲嫩草| 欧美性长视频在线观看| 69av精品久久久久久| 国产乱人伦免费视频| 亚洲国产精品成人综合色| 亚洲va日本ⅴa欧美va伊人久久| 两个人视频免费观看高清| 成人18禁在线播放| 免费av毛片视频| 黄色女人牲交| 国产成人欧美在线观看| 此物有八面人人有两片| 99riav亚洲国产免费| 777久久人妻少妇嫩草av网站| 亚洲aⅴ乱码一区二区在线播放 | 老汉色av国产亚洲站长工具| 久久青草综合色| 啦啦啦 在线观看视频| 麻豆一二三区av精品| 久久久国产精品麻豆| 韩国av一区二区三区四区| 国产伦人伦偷精品视频| 欧美激情久久久久久爽电影 | 欧美精品啪啪一区二区三区| 中文字幕另类日韩欧美亚洲嫩草| 不卡av一区二区三区| 人成视频在线观看免费观看| 久热爱精品视频在线9| 在线观看免费日韩欧美大片| 黄色 视频免费看| 别揉我奶头~嗯~啊~动态视频| 女人高潮潮喷娇喘18禁视频| xxx96com| 欧美成人性av电影在线观看| 女人被躁到高潮嗷嗷叫费观| 法律面前人人平等表现在哪些方面| 亚洲精品粉嫩美女一区| 激情在线观看视频在线高清| 人人妻人人澡人人看| av视频免费观看在线观看| 一级毛片高清免费大全| 一边摸一边抽搐一进一小说| 91成人精品电影| 国产成人av激情在线播放| 狠狠狠狠99中文字幕| ponron亚洲| 成人18禁在线播放| 国产精品精品国产色婷婷| 国产视频一区二区在线看| 国产97色在线日韩免费| www.999成人在线观看| 香蕉久久夜色| 电影成人av| 国产视频一区二区在线看| 麻豆国产av国片精品| 亚洲精品在线观看二区| 久久久久久亚洲精品国产蜜桃av| 亚洲成av片中文字幕在线观看| 欧美激情久久久久久爽电影 | 日韩 欧美 亚洲 中文字幕| 免费一级毛片在线播放高清视频 | 97碰自拍视频| 最近最新中文字幕大全电影3 | 高清毛片免费观看视频网站| 亚洲国产精品sss在线观看| 国内精品久久久久精免费| 啪啪无遮挡十八禁网站| 男女做爰动态图高潮gif福利片 | 国产一卡二卡三卡精品| 午夜福利影视在线免费观看| 桃红色精品国产亚洲av| 国产精品久久久久久亚洲av鲁大| 成人亚洲精品av一区二区| 88av欧美| 亚洲一区二区三区不卡视频| 老鸭窝网址在线观看| 神马国产精品三级电影在线观看 | www.精华液| or卡值多少钱| 欧美日韩黄片免| 色综合欧美亚洲国产小说| 成人亚洲精品一区在线观看| 欧美成人午夜精品| 亚洲成国产人片在线观看| 国产黄a三级三级三级人| 又紧又爽又黄一区二区| 可以免费在线观看a视频的电影网站| 久久久水蜜桃国产精品网| 久久国产精品影院| 老司机午夜福利在线观看视频| 国产亚洲欧美在线一区二区| 国语自产精品视频在线第100页| 国产区一区二久久| 19禁男女啪啪无遮挡网站| 一边摸一边做爽爽视频免费| 亚洲欧美精品综合久久99| 一个人免费在线观看的高清视频| 亚洲中文字幕一区二区三区有码在线看 | 亚洲三区欧美一区|