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

    應(yīng)對水文序列非一致性變化影響的溯源重構(gòu)法研究

    2021-08-20 07:13:34毅,李
    水利學報 2021年7期
    關(guān)鍵詞:洪峰流量水文一致性

    秦 毅,李 時

    (省部共建西北旱區(qū)生態(tài)水利國家重點實驗室西安理工大學,陜西西安 710048)

    1 研究背景

    近年來,大量研究成果表明實測水文序列存在明顯的非一致性特征[1-4]。它帶來兩方面問題:(1)在非一致性條件下,繼續(xù)使用以一致性樣本為基礎(chǔ)的預(yù)測分析方法,包括傳統(tǒng)水文頻率分析方法、多元統(tǒng)計分析的Coupla方法、機器學習等方法是否能獲得合理的成果,以及對未來的預(yù)測應(yīng)該如何評價。(2)水文模型參數(shù)的率定工作變得困難,因為參數(shù)率定采用資料的非一致性導致了不確定性的增加。相比較而言,第一個問題在設(shè)計分析方面顯得更為突出,例如水利工程、水土保持工程及生態(tài)工程的設(shè)計,洪旱災(zāi)害的防治等都離不開對實測樣本序列的統(tǒng)計分析與預(yù)測,除非研究對象的物理機制被真實地掌握,否則所做出的水文分析與預(yù)測的結(jié)果只有在一致性(平穩(wěn)性)序列條件下才有意義。因此學者們投入了大量研究,包括:(1)對水文序列非一致性變異診斷的深入研究,給出了相對成熟的診斷方法:如相關(guān)系數(shù)檢驗法、Spearman秩次相關(guān)法、Kendall 秩次相關(guān)法、Mann-Kendall檢驗法[5-6]、R/S分析方法[7]、非參數(shù)Pettitt檢驗法[8]等。在此基礎(chǔ)上,謝平等[9]構(gòu)建了水文變異診斷系統(tǒng),更全面地反映了水文序列的變異特征。(2)解決非一致性條件下的頻率分析問題,這是工程設(shè)計迫切需要解決的問題[10]。為應(yīng)對非一致性序列的影響,除了采用水文分析計算相關(guān)規(guī)范提出的還原或還現(xiàn)方法外,近年來有學者們將研究焦點放在了時變矩法[11-12]上,并據(jù)此建立時變概率分布,再通過重現(xiàn)期或可靠度與時變概率分布間的關(guān)系確定出工程設(shè)計值。當前比較有代表性的方法包括基于修正重現(xiàn)期的EWT[13-14]和ENE方法[15-16],和基于可靠度的ER[17]、DLL[18]和ADLL[19]等方法。

    然而,這些方法仍然存在一定的爭議。首先是對趨勢的看法,有學者認為觀測到的水文序列的趨勢并非是非一致性特征的表現(xiàn),而可能是平穩(wěn)過程中的周期性頻率波動所致[20],甚至認為趨勢這個詞都缺少明確的定義[21];再者,實踐中,在工程設(shè)計壽命內(nèi),一定可靠度下的設(shè)計值將隨著計算時段的起始時間以及統(tǒng)計參數(shù)與協(xié)變量間關(guān)系擬合曲線的類型而變化,這就意味著未來設(shè)計值的可靠性嚴重依賴統(tǒng)計參數(shù)的時變特征,而這種時變特征卻因為時間序列遍歷性的缺失,造成對未來水文序列統(tǒng)計參數(shù)預(yù)測的不確定性[22]。正如Serinaldi等[23]指出的,如果模型結(jié)構(gòu)無法通過演繹推斷的方式獲得,特別是描述分布參數(shù)非平穩(wěn)變化模型是僅通過歸納推理擬合得來,則一個額外不確定性來源被無形中引入了模型結(jié)構(gòu)中,使得所得到的非平穩(wěn)模型無法在實際應(yīng)用中提高設(shè)計值的可信度與準確性,而模型的這種不確定性將很容易導致不具物理意義的結(jié)果產(chǎn)生;此外,利用時變參數(shù)識別非平穩(wěn)序列的分布函數(shù)過程中,所采用的經(jīng)驗頻率計算方法也是引起不確定性的重要因素,因為計算經(jīng)驗頻率的方法也同樣要求樣本為簡單隨機樣本,基于非平穩(wěn)序列直接進行經(jīng)驗頻率計算顯然不恰當??梢?,對于統(tǒng)計特征隨時間變化的單一隨機變量,其分布研究尚存在如此多的問題,更何況解決多元時變隨機變量的聯(lián)合分布問題,難度可想而知。

    很顯然,造成“難度”的核心是樣本序列的非簡單性。如果能夠?qū)⒎且恢滦詷颖拘蛄修D(zhuǎn)換為一致性序列,則因為對此簡單序列已經(jīng)有了成熟的分析理論、技術(shù)和用好技術(shù)的方法,前述困難將不復(fù)存在。因此,本文提出溯源重構(gòu)法,它是基于反映水文現(xiàn)象物理機理的源函數(shù)將非一致性序列重構(gòu)為新的一致性序列的方法。依靠重構(gòu)的一致性序列,解決需要基于一致性序列才能進行分析工作的問題。事實上,當前已存在將非一致性序列轉(zhuǎn)變?yōu)橐恢滦孕蛄械姆椒?,即還原或還現(xiàn)法、分解合成法,本文將在后面就這兩種方法與溯源重構(gòu)法進行比較與討論。

    2 溯源重構(gòu)法的理論與方法

    2.1 源函數(shù)的概念一般來說,一個水文變量Y的變化是多個影響因子X1,…,XN共同作用的結(jié)果,見圖1。在排除其他影響因子影響的條件下,單獨考察影響因子Xi對水文變量Y的作用,將這種作用規(guī)律記為fi(Xi),則fi(Xi)是Xi對Y的物理作用機制,因此這里稱fi(Xi)是Xi作用于Y的源函數(shù)。作為物理機制,源函數(shù)是不會改變的,就像河道流量Q=AV,其中A是過流斷面的面積,V是流速,則V和A各自對Q的作用規(guī)律或說機制,也即源函數(shù)就是fV(V)=V,fA(A)=A。當V或A或兩者共同隨時間變化時,流量也會隨著V(t)或A(t)或兩者在不同時刻的取值而表現(xiàn)出隨時間變化的Q(t)。但是,不論V和A怎樣取值,源函數(shù)fV(V)=V,fA(A)=A保持不變,而流量等于這兩個源函數(shù)相乘的機制也不會改變。即物理機制不隨時間變化,變化的永遠是Xi的狀態(tài)。

    圖1 各影響因子對研究變量Y的源函數(shù)作用關(guān)系

    嚴格的說,源函數(shù)的確定有兩種方法:理論推導(數(shù)學解析)和實驗分析。實際上,由于水文現(xiàn)象變化的復(fù)雜性,即便我們認可源函數(shù)的存在,也會因為人類對自然的有限認知而很難得到絕對準確的物理機制數(shù)學表達,至少當前很少看到對單一因子影響的研究。由于事物之間的作用規(guī)律可以隱含于統(tǒng)計規(guī)律中,所以當下可以借用Y和Xi間的統(tǒng)計關(guān)系來近似估計源函數(shù)。

    既然源函數(shù)或由源函數(shù)構(gòu)成的函數(shù)不隨時間變化,我們就可以利用這個特征將非一致性序列轉(zhuǎn)變?yōu)橐恢滦孕蛄小?/p>

    2.2 溯源重構(gòu)方法的基礎(chǔ)模型欲利用源函數(shù)時不變的特點,就要表達出多個影響因子各自對Y作用后的一般綜合關(guān)系。鑒于數(shù)學解析方法暫時行不通,故采用建立統(tǒng)計模型的方法來表征研究對象Y與解釋變量源函數(shù)之間的關(guān)系。眾所周知,數(shù)學的基本運算包括加(減)與乘(除),結(jié)合運算法則,Y與解釋變量源函數(shù)之間的關(guān)系有兩種一般性廣義表達方式,一是各解釋變量間相互獨立時采用的疊加模式:

    另一個是可以考慮解釋變量間存在相關(guān)性的乘積模式:

    其中:N指作用于研究變量Y的所有影響因子的個數(shù),這些影響因子包括人類認識到的和尚未認識的;fi(Xi(t))為影響因子Xi的源函數(shù)。假設(shè)在N個影響因子中,有n個因子是人類目前已認識到的,它們的“源”函數(shù)分別為f1(X1(t)),…,fn(Xn(t));未知影響因子及其“源”函數(shù)分別為fn+1(Xn+1(t)),…,fN(XN(t))。則式(1)和式(2)改寫為:

    分別記δ=fn+1(Xn+1(t))+…+fN(XN(t))和δ*=fn+1(Xn+1(t))…fN(XN(t)),他們代表人類未知的影響因素對Y的作用,通??梢哉J為是自然隨機變量,具有概率分布,并具有對應(yīng)的均值μ、μ*,和方差σ2、σ*2。令Ks(t)=f1(X1(t))+…+fn(Xn(t))和Kp(t)=f1(X1(t))…fn(Xn(t))。在生產(chǎn)實踐當中,由于我們更關(guān)心影響因子取給定值時,Y的變化情況,與回歸分析的考慮相同,在計算Y之前,解釋變量X1(t),…,Xn(t)的取值已經(jīng)確定,可以作為可控變量處理[24]。于是Ks(t)和Kp(t)在t時刻為常數(shù),故Ks(t)和Kp(t)是已知的時間函數(shù)。重新寫上述式(3)、式(4)為:

    對比兩種表達式下Y(t)的條件數(shù)學期望和條件方差如下。

    疊加模型:

    乘積模型:

    不難看到,兩種模型下Y(t)的條件均值均是時變函數(shù),一旦Xi發(fā)生趨勢性變化,Y(t)的條件數(shù)學期望也就出現(xiàn)趨勢性變化;與均值不同,在疊加模型下的方差為常數(shù),而乘積模型下方差是時變函數(shù),且會隨Xi的趨勢性變化而產(chǎn)生更強的趨勢性。

    現(xiàn)在來考察水文現(xiàn)象變化的普遍特點。水文現(xiàn)象是在水文系統(tǒng)內(nèi)各要素間普遍存在相互影響的背景下發(fā)生的。各解釋變量間基本不存在絕對的獨立性,他們的相互作用使水文系統(tǒng)呈現(xiàn)出高度復(fù)雜的自適應(yīng)非線性狀態(tài)。隨著近年來的環(huán)境變化,大量的研究成果表明不僅水文變量的量級表現(xiàn)出了非平穩(wěn)變化,而且水文極端事件也急劇增多[25-26]。從統(tǒng)計的角度來看,這種現(xiàn)象表明了水文序列的一、二階矩均發(fā)生了變化[27-28],特別是極端事件的極端特性越強,二階矩的變化就越不能被忽視[29]。因此,可以得到結(jié)論,乘積模式相較于疊加模式更適用于描述具有復(fù)雜相關(guān)性的水文系統(tǒng)。故這里采用乘積模型(11)作為溯源重構(gòu)法的基礎(chǔ)模型。

    2.3 溯源重構(gòu)函數(shù)與重構(gòu)序列考慮到Y(jié)(t)的非一致性變化是因為解釋變量的非一致性變化引起,那么如果將所有帶有非一致性變化的影響因子的作用從Y(t)中剔除,見式(12),則剩余量將表現(xiàn)出一致性。

    在式(12)中,m是非一致性變化的影響因子個數(shù)。其右端為Y(t)剔除全部非一致性因子影響后剩余一致性影響因子以及尚未被人類認知的影響因子源函數(shù)的乘積,它是平穩(wěn)的隨機過程。如果只剔除了部分非一致性變化因子X1,X2,…,Xl(l

    理論上,當引起Y非一致性變化的影響因子全部被識別出來,即l=m,且這些因子的源函數(shù)fi(Xi(t))已有準確的表達,那么由溯源重構(gòu)函數(shù)計算出來的溯源重構(gòu)序列RSt={rs1,rs2,…,rsn}為一致性序列,可用于需要一致性序列才能進行工作的情況。生產(chǎn)實踐中,因人類認知和現(xiàn)有方法的局限,我們不可能得到絕對的平穩(wěn)過程,這也不是我們的目標,但隨著非一致性變化的解釋變量對應(yīng)的正確源函數(shù)不斷加入溯源重構(gòu)函數(shù)中,則溯源重構(gòu)序列將逐漸趨向并最終趨向于平穩(wěn),接近于一個隨機噪聲。

    需要提醒的是,由于估計源函數(shù)存在不確定性,所以使用中必須配合溯源重構(gòu)序列RSt的一致性檢驗。若檢驗結(jié)果為一致,則同時接受所確定的影響因子及其源函數(shù)的估計形式,否則,仍需進一步研究影響因子及其源函數(shù)。

    由于研究變量Y和影響因子間的因果關(guān)系,研究變量Y的各種非一致性變化(線性、非線性、突變)與影響因子的相應(yīng)非一致性變化一致,這使得溯源重構(gòu)法去除非一致性時并不受非一致性變化類型的影響,即具有普遍性。相較于對非一致性變化的Y進行純數(shù)學模擬的時變矩法,溯源重構(gòu)方法更注重找到引起Y變化的非一致性變化影響因子Xi和它對Y的作用規(guī)律,并以此規(guī)律將非一致性變化從Y序列中排除,這正是溯源重構(gòu)法能夠?qū)⒎且恢滦宰兓腨(t)變?yōu)橐恢滦孕蛄蠷S(t)的關(guān)鍵所在。

    3 實例研究

    3.1 研究區(qū)域與數(shù)據(jù)利用陜西省榆林市佳蘆河流域申家彎水文站年最大洪峰系列資料作為實例研究,考察溯源重構(gòu)法構(gòu)建一致性序列的效果。由于溯源重構(gòu)法的重點是依據(jù)非一致性變化影響因子對水文變量的作用(即源函數(shù))來確定溯源重構(gòu)函數(shù)的,因此在研究中,首先應(yīng)研究確定非一致性變化的影響因子。實踐中,可以依據(jù)水文現(xiàn)象的相關(guān)理論知識并結(jié)合實地踏勘調(diào)研等手段,初步確定研究變量的主要影響因子,之后對這些主要影響因子進行趨勢分析,找出造成水文變量非一致性變化的影響因子。

    佳蘆河位于黃河中游的河龍區(qū)間,毛烏素沙漠南緣,陜北黃土高原北段,屬于黃土高原溝壑區(qū)。河流由西北流向東南,至佳縣佳蘆鎮(zhèn)木廠灣村注入黃河。河流全長93 km,流域面積1134 km2,出口控制站為申家灣水文站。流域洪水由降雨引發(fā)。由于氣候較為干旱,自然條件較差,該地區(qū)容易發(fā)生風沙災(zāi)害,水土流失現(xiàn)象普遍。為了控制水土流失,流域內(nèi)修建了數(shù)十座淤地壩。淤地壩一般由兩大件組成——壩和泄洪臥管。為了使入庫泥沙得到沉積,實現(xiàn)多攔沙少攔水的功效,淤地壩運行方式設(shè)計成:發(fā)生洪水時,高含沙洪水首先被攔截在庫內(nèi),后由臥管分數(shù)天將洪水排出庫外,且要求臥管排水速度緩慢,排水流量一般為1~1.5 m3/s。結(jié)果導致出口斷面洪峰較建壩前明顯減小,庫壩泄流對流域出口的洪峰流量的貢獻可以忽略不計。這種效果等同于流域產(chǎn)生洪峰的產(chǎn)流面積減少,減少的數(shù)量為庫壩控制面積。這里稱流域面積與庫壩控制面積的差為有效產(chǎn)流面積;另一方面退耕還林以及淤地壩的蓄水,使得坡面較大范圍內(nèi)的土壤含水率增加,促進了植被增長,導致流域植被截留、下滲、匯流等條件均發(fā)生改變,也造成了洪峰流量序列的非一致性變化。因此,選擇有效產(chǎn)流面積與植被覆蓋度作為佳蘆河流域年最大洪峰流量序列非一致性變化的主要影響因子。

    洪峰流量序列采用申家灣水文站1957—2010年實測數(shù)據(jù)進行分析。1957—2010年有效產(chǎn)流面積序列則根據(jù)流域面積與逐年修建淤地壩的控制面積計算確定。植被覆蓋度(FVC)是利用landsat4-5 TM 30 m×30 m遙感資料,通過歸一化植被指數(shù)(NDVI)的分析計算得到[30]。由于植被覆蓋度序列長度受限于遙感資料的年限,僅有1988—2010年23年資料。但在1980年代前,流域人類活動影響較弱,且退耕還林等生態(tài)工程尚未起步,因此1980年代以前的植被狀態(tài)較為穩(wěn)定,本文將1988年前的植被覆蓋度均取值為1988年數(shù)據(jù)。

    圖2和圖3分別展示了佳蘆河流域出口斷面水文站申家灣站的年最大洪峰流量序列的變化和流域有效產(chǎn)流面積的變化??梢钥吹剑饔蚝榉辶髁亢陀行Мa(chǎn)流面積在1975年左右均發(fā)生了突變。圖4顯示了流域植被覆蓋度的變化情況。對各序列采用Pearson 相關(guān)系數(shù)檢驗判斷序列的線性趨勢,Spearman 秩相關(guān)系數(shù)檢驗判斷包括非線性趨勢在內(nèi)的一階矩趨勢,采用Breusch-Pagan 檢驗[27]判斷二階矩的趨勢。因1988年前的植被覆蓋度數(shù)據(jù)為作者外延處理確定的,對植被覆蓋度趨勢性檢驗時采用1988—2010年的數(shù)據(jù)進行分析。由于淤地壩因逐年修建而增多,且建壩本身不屬于隨機事件,所以有效產(chǎn)流面積序列呈現(xiàn)顯著單調(diào)下降的趨勢,也沒有一般隨機序列的波動變化。

    圖3 佳蘆河流域有效產(chǎn)流面積序列

    圖4 佳蘆河流域植被覆蓋度序列

    檢驗結(jié)果(表1)表明年最大洪峰流量序列的均值和方差都出現(xiàn)了顯著減少趨勢,但植被覆蓋度的均值出現(xiàn)顯著上升趨勢。由于年最大洪峰流量一般在汛期發(fā)生,本文對佳蘆河流域申家灣站的各時段致洪降雨序列進行了趨勢性檢驗??紤]到黃土丘陵溝壑區(qū)降雨過程一般歷時較短,不超過12 h,因此最大12 h降雨能夠代表相應(yīng)洪水的致洪暴雨,故在此處僅展示最大12 h降雨變化情況(圖2)。同時將與佳蘆河流域較近的榆林氣象站最大日降雨序列作為參考,一并進行檢驗。檢驗結(jié)果(表2)顯示,佳蘆河流域致洪降雨序列未發(fā)生顯著變化(包括線性和非線性一階矩、二階矩)。因此排除了降雨變化造成洪峰流量非一致性變化的可能。

    圖2 申家灣站年最大洪峰流量序列與致洪暴雨序列

    表1 佳蘆河流域洪峰流量序列及其影響因子的趨勢性分析

    表2 佳蘆河流域致洪降雨序列趨勢性分析

    3.2 溯源重構(gòu)序列的生成通過對流域?qū)崪y資料的分析,我們可以得到這樣的認識:佳蘆河流域的降雨未發(fā)生顯著的非一致性變化,申家灣站年最大洪峰流量序列的非一致性變化主要產(chǎn)生于流域有效產(chǎn)流面積和植被覆蓋度的顯著變化。因此,根據(jù)溯源重構(gòu)法的思想,重構(gòu)時不考慮降雨因子,僅考慮有效產(chǎn)流面積與植被覆蓋度兩個因子。在對洪峰流量序列進行重構(gòu)之前,首先要分析出有效產(chǎn)流面積和植被覆蓋度對洪峰流量作用的源函數(shù)。

    3.2.1 源函數(shù)的建立方法 源函數(shù)是溯源重構(gòu)法的關(guān)鍵,也是提升對水文過程認知的重要基礎(chǔ),是不應(yīng)該被忽略的學術(shù)問題。它的確定應(yīng)該來自于理論推導或?qū)嶒?,但現(xiàn)階段,介于水文過程的復(fù)雜性,基于理論和實驗手段成功確定出的源函數(shù)幾乎看不到。因此,本著影響因子的作用規(guī)律必然體現(xiàn)在統(tǒng)計規(guī)律中的認識,我們可以通過回歸的方法來逐個估計源函數(shù)的近似值。為敘述方便,這里令Y表示洪峰流量,X1和X2分別表示流域有效產(chǎn)流面積和植被覆蓋度。源函數(shù)建立的具體做法:首先建立第一個因子X1與Y的關(guān)系f1(X1(t)),并作為X1的源函數(shù)作用于Y,之后將該因子的影響從Y序列中排除,即,再建立剩余序列與第二個因子X2的源函數(shù)f2(X2(t))。推廣而言,若選擇了m個有趨勢變化的影響因子,則第3個影響因子的源函數(shù)由與第3個因子的關(guān)系估計,依次進行,直到m個因子的源函數(shù)被估計出。

    顯然,這種做法會帶來因影響因子的引入順序改變各因子源函數(shù)形式的問題。理論上,真實的源函數(shù)將不會受重構(gòu)順序的干擾。而實際中,由于源函數(shù)是估計函數(shù),重構(gòu)順序的干擾是可能存在的,為避免這個問題,可根據(jù)各影響因子與研究變量之間相關(guān)性的強弱順序確定重構(gòu)順序。

    3.2.2 有效產(chǎn)流面積與植被覆蓋度的單因子源函數(shù) 所謂單因子指只有一個非一致性性影響因子。眾多的研究成果證明洪峰流量和流域面積之間呈現(xiàn)冪函數(shù)關(guān)系[31-32]。權(quán)琦澤[33]的研究結(jié)果表明榆林地區(qū)洪峰流量與流域面積的0.73次方成正比,故假定有效產(chǎn)流面積的源函數(shù)的估計為

    植被對洪峰流量的作用關(guān)系研究不多見。根據(jù)周宏飛等[34]在徑流試驗場的實驗研究成果,徑流深與植被覆蓋度FVC間呈現(xiàn)指數(shù)函數(shù)關(guān)系:

    式中:R為徑流深,mm;k為經(jīng)驗參數(shù)。借助洪峰與洪量間的冪函數(shù)關(guān)系、洪量與徑流深間的關(guān)系,便可由式(15)得到洪峰流量與植被覆蓋度間的作用規(guī)律:Qm∝eαFVC,因此植被覆蓋度FVC對洪峰流量作用的源函數(shù)估計式為式(16):

    通過實測洪峰流量和eαFVC之間的回歸,得佳蘆河流域的參數(shù)α為-0.071,即:

    3.2.3 年最大洪峰流量的溯源重構(gòu)序列 當確立了有效產(chǎn)流面積和植被覆蓋度這兩個變化因子的源函數(shù)后,根據(jù)式(13)定義的溯源重構(gòu)函數(shù),建立佳蘆河流域的溯源重構(gòu)函數(shù)。為了比較溯源重構(gòu)方法的有效性,分別建立單個因子和兩因子的溯源重構(gòu)函數(shù)如下:

    利用有效產(chǎn)流面積的單因子重構(gòu):

    利用植被覆蓋度的單因子重構(gòu):

    利用有效產(chǎn)流面積和植被覆蓋度兩因子重構(gòu)(方法見3.2.1節(jié)):

    其中fFVC(FVC)=eαFVC中的α=-0.059則由RSA,t與eαFVC的回歸得到。

    理論上講,若源函數(shù)是通過理論分析與實驗分析的方式確定,則無論是單因子重構(gòu)還是兩因子重構(gòu),源函數(shù)不應(yīng)發(fā)生變化。這里α 取值的差異正是源函數(shù)的估計受到重構(gòu)順序影響的體現(xiàn)。因此,實際應(yīng)用中建議根據(jù)各影響因子與研究變量之間相關(guān)性的強弱順序確定重構(gòu)順序。將洪峰流量序列Qm,t及其影響因子序列Ae,t和FVCt的值分別代入式(18)—(20)中,得到相應(yīng)的溯源重構(gòu)序列RSA,t、RSFVC,t和RSdouble,t。

    4 討論

    按照本文提及的溯源重構(gòu)理論建立的溯源重構(gòu)方法,其實質(zhì)是一種轉(zhuǎn)換方法,所獲得的重構(gòu)序列RSt理應(yīng)是一個具有一致性的隨機序列。作為驗證,本文以第3節(jié)中的溯源重構(gòu)序列RSA,t、RSFVC,t和RSdouble,t為對象,通過考察其趨勢性變化,探究溯源重構(gòu)理論的正確性和溯源重構(gòu)序列的合理性。同時,鑒于以下兩個原因,對重構(gòu)序列的特征進行探討同樣顯得必要。①雖然RSt(如RSA,t、RSFVC,t和RSdouble,t)不是直接觀測變量,但因它是由觀測變量計算出來,因此有物理意義:相當于模數(shù),應(yīng)該具有它的特征;②實踐中,在溯源重構(gòu)時,隨趨勢性影響因子識別的完整性不同,所獲得的重構(gòu)序列會具有不同的特征。因此重構(gòu)序列的特征分析對于溯源重構(gòu)序列RSt的應(yīng)用起著理論支撐作用。

    4.1 溯源重構(gòu)序列的趨勢特征采用同樣的假設(shè)檢驗方法,對上述溯源重構(gòu)序列的趨勢性加以檢驗,檢驗結(jié)果見表3。可以看出,具有顯著趨勢性變化的原洪峰流量序列經(jīng)單因子溯源重構(gòu)后,重構(gòu)序列的各種相關(guān)系數(shù)均比原序列相應(yīng)相關(guān)系數(shù)減小,說明不論是線性趨勢還是非線性趨勢均得到一定的消除;特別是雙因子重構(gòu)后,Pearson相關(guān)系數(shù)和Spearman秩相關(guān)系數(shù)均小于相應(yīng)的臨界值,說明溯源重構(gòu)序列RSdouble,t已無顯著趨勢存在,且Breusch-Pagan 檢驗統(tǒng)計量c也由25.406 減小到5.157,盡管仍大于臨界值,但洪峰流量序列二階矩的非平穩(wěn)性得到很大改善,說明引入非平穩(wěn)變化的影響因子越合理,重構(gòu)序列的平穩(wěn)性越好。不僅如此,原序列存在的突變也得到消除,見圖5。

    圖5 溯源重構(gòu)序列與年最大洪峰流量序列的對比

    表3 年最大洪峰流量溯源重構(gòu)序列一致性檢驗

    檢驗結(jié)果表明,溯源重構(gòu)方法在消除非一致性序列的非一致性(包括線性、非線性趨勢和突變)方面具有統(tǒng)計意義上的有效性。

    4.2 溯源重構(gòu)序列的分布特征分布特征是隨機變量取值特性的完整描述,是研究隨機現(xiàn)象的核心。這里討論年最大洪峰流量溯源重構(gòu)序列的分布函數(shù),并與基于時變矩獲得的洪峰流量序列(原序列)的分布函數(shù)進行比較。同時,針對當前工程師們無奈之選,即無視非平穩(wěn)的存在,直接進行傳統(tǒng)頻率分析得出的分布函數(shù)也一并列入比較。

    針對3.2 節(jié)得到的佳蘆河年最大洪峰流量的溯源重構(gòu)序列RSdouble,t,本文選取水文領(lǐng)域常用的Weibull(WEI),Log-normal(LN),Gamma(GA),Gumbel(GU)和GEV 五種分布類型作為備選分布,利用線性矩法[35]估計分布參數(shù),依據(jù)納什效率系數(shù)[36-37]、均方根誤差RMSE等擬合優(yōu)度準則確定最優(yōu)分布。申家灣站溯源重構(gòu)序列的最優(yōu)分布為GEV,對應(yīng)的分布參數(shù)分別為μ=9.657,σ=83.525,κ=0.878。對于非一致性序列的時變分布則基于GAMLSS模型[38]采用時變矩法確定,備選分布仍選取上述五種類型,以有效產(chǎn)流面積與植被覆蓋度作為協(xié)變量,評價準則為AIC準則[39],具體見表4。而無視非一致性存在,直接以非一致性原洪峰流量序列估計的總體分布也列于表4。

    表4 不同方法確定的年最大洪峰流量序列的分布類型

    理論上,年最大洪峰流量序列屬于極值序列,它的分布應(yīng)該是極值型。從表4可以看到不同屬性序列和不同分析方法所挑選出的分布函數(shù)不同,其中只有溯源重構(gòu)序列的分布符合極值型,本文作者對陜西省榆林地區(qū)大理河等6個流域的年最大洪峰流量溯源重構(gòu)序列進行分布函數(shù)分析,結(jié)果均為極值類分布,這說明溯源重構(gòu)法可以體現(xiàn)原序列的極值屬性。這種既平穩(wěn)又符合抽樣特點的序列表明溯源重構(gòu)法的理論正確,所得序列可用。對于非一致性原序列,相較于傳統(tǒng)方法,時變矩法因考慮了原序列非一致性的特征,所以挑選出了與極值類型相接近的GA分布。有趣的是若直接應(yīng)用非一致性原序列,采用傳統(tǒng)方法分析出的分布函數(shù)卻是不同于極值型分布的LN分布。這一結(jié)果從側(cè)面反映了水文序列分布的確定會受到參數(shù)估計方法的影響,即估計方法越適應(yīng)序列屬性,所選擇出的分布越可靠。

    利用頻率分析推求水文設(shè)計值是水資源利用、管理與保護領(lǐng)域的最基本問題。合理的設(shè)計值理應(yīng)由合理的分布確定。Li等[22]在研究淤地壩影響下的設(shè)計洪水計算問題時,所采用的去非平穩(wěn)性方法實質(zhì)上正是溯源重構(gòu)思想的體現(xiàn)。其獲得的合理結(jié)果同樣證明溯源重構(gòu)方法獲得的概率分布合理,同時說明該方法在非一致性序列頻率分析方面具有良好的適應(yīng)性。

    4.3 與其他構(gòu)建一致性序列方法比較溯源重構(gòu)法的思想是通過溯源重構(gòu)將非一致性序列轉(zhuǎn)變?yōu)橐恢滦孕蛄?,為以一致性序列作為基礎(chǔ)的分析計算工作提供樣本資料。這種想法同樣體現(xiàn)在還原或還現(xiàn)法和分解合成法中。但由于方法所依據(jù)的原理不同,方法的應(yīng)用效果隨之不同。

    還原或還現(xiàn)法依據(jù)的原理是通過對水文變量系統(tǒng)性變化的增量回加到原序列或?qū)λ淖兞窟M行修正,達到使水文變量序列一致的目的。在增量的量化或修正值的量化工作中,由于要求資料的種類多、數(shù)量大,且對資料的準確性依賴程度高,使用不易。同時牽扯的一些概念難以精確界定[40],存在方法失真或失效[41]等問題。即便用分布式水文模型直接計算增量,也會面臨資料的非一致性導致模型參數(shù)率定(優(yōu)化法)的不確定性加大問題[41]。正因如此,本文不做具體的分析對比。

    分解合成法[42]認為,“無論水文現(xiàn)象的變化多么復(fù)雜,水文序列總可以分解成確定性成分和隨機性成分這兩部分。一般來說,水文序列的確定性成分主要受人類活動的影響,但并不排除氣候因素(如氣候轉(zhuǎn)型期)和下墊面因素(如火山爆發(fā)、地震等)的影響,其變化規(guī)律可以在較短的工程年代里發(fā)生緩慢的漸變或劇烈的突變,因此水文序列中確定性成分的變化規(guī)律往往是非一致的。而水文序列的隨機性成分則主要受氣候、地質(zhì)等因素的影響,其變化規(guī)律需要一個漫長的地質(zhì)年代才能改變,因此水文序列中隨機性成分的統(tǒng)計規(guī)律是相對一致的。”于是謝平等[42]建議采用成因分析法與統(tǒng)計分析法分別對確定性成分和隨機性成分進行識別與檢驗,并對確定性成分進行擬合計算,對隨機性成分采用傳統(tǒng)方法進行頻率分析計算。后者意味著分解合成法獲得的隨機變量樣本是一個簡單隨機樣本。下面將繼續(xù)以佳蘆河流域年最大洪峰流量序列資料,研究分析分解合成法生成的隨機樣本特征。

    根據(jù)分解合成法原理,非一致性水文序列Xt可按下式表示:

    式中:Yt為確定的非周期成分(包括趨勢、跳躍等);Pt為確定性的周期成分(包括簡單的或復(fù)合周期的成分等);St為隨機成分。

    為避免周期對趨勢的影響,這里首先進行周期分析,利用小波分析[43]得到周期序列存在21年、9年兩個顯著周期,由此獲得周期函數(shù)如式(22):

    剔除周期后發(fā)現(xiàn)序列仍存在明顯的趨勢,通過回歸分析得到式(23)所示的趨勢成分:

    根據(jù)式(21)剔除周期成分與趨勢成分后得到隨機序列St,見圖6,對該隨機序列進行趨勢檢驗,結(jié)果見表5。

    圖6 分解合成法構(gòu)建的佳蘆河年最大洪峰流量一致性序列

    表5 分解合成法構(gòu)建的佳蘆河年最大洪峰流量序列的一致性檢驗

    分解合成法獲得的佳蘆河年最大洪峰流量的重構(gòu)序列十分有效地消除了一階矩的顯著線性趨勢,但Spearman秩相關(guān)系數(shù)依然大于臨界值,表明其一階矩的非線性型非一致性仍然顯著存在;方差的非一致性檢驗統(tǒng)計量χ為23.064,幾乎沒有任何改善,說明分解合成法重構(gòu)的隨機成分并非簡單樣本。

    產(chǎn)生上述現(xiàn)象的原因可以概化為兩個方面。一是由于定義兩種成分的物理機制均不夠清晰,例如分解合成法原理述及的“水文序列的確定性成分主要受人類活動的影響,但并不排除氣候因素(如氣候轉(zhuǎn)型期)和下墊面因素(如火山爆發(fā)、地震等)的影響”與“隨機性成分則主要受氣候、地質(zhì)等因素的影響,其變化規(guī)律需要一個漫長的地質(zhì)年代才能改變”,其中氣候轉(zhuǎn)型期以及需要漫長的地質(zhì)年代才能改變的氣候因素無法清晰界定,更不能明確表達,因此確定性成分的擬合具有較大的不確定性,由實測資料推求出的隨機性成分為實測資料與確定性成分之差,就會具有較大的不確定性;二是該方法假定水文序列Xt的各個成分滿足線性疊加特性(即疊加模型)[44]。因此,根據(jù)上文疊加模型與乘積模型的對比分析,分解合成法構(gòu)建的“相對一致”性序列僅能實現(xiàn)均值一致,而方差非一致的現(xiàn)象就不難理解了。這也成為了疊加模型無法描述序列的二階矩變化的例證。

    相比之下,溯源重構(gòu)法認為水文變量Y的變化是一系列影響因子自身的非一致性變化引起的。這些因子的非一致性影響通過它們各自的源函數(shù)作用于Y。如果將全部非一致性變化因子的作用從Y序列中剔除,則剔除后的序列RSt就是一致性序列。與分解合成法不同,溯源重構(gòu)法中的隨機成分是指未被人類認知的次要因子的作用,無需專門率定。同時溯源重構(gòu)法的基礎(chǔ)模型采用乘法模式,更符合水文系統(tǒng)的非線性相依特征。溯源重構(gòu)法的物理含義明確,分析的不確定性相對較小。

    5 結(jié)論

    眾所周知,任何現(xiàn)象都有其因果關(guān)系,以數(shù)學術(shù)語來說,即因變量與自變量的關(guān)系。作者認為,水文變量Y的非一致性變化是由于自變量(影響因子)自身的非一致性變化所引起的。據(jù)此提出溯源重構(gòu)法,即從源函數(shù)概念出發(fā),依據(jù)因變量與自變量間的物理機理所建立起的統(tǒng)計模型,從因變量Y序列中剔除非一致性變化的自變量影響,實現(xiàn)非一致性序列Yt向新的一致性序列RSt的轉(zhuǎn)化,從而可以服務(wù)于以一致性序列為基礎(chǔ)的分析計算工作。溯源重構(gòu)法在統(tǒng)計原理的基礎(chǔ)上充分考慮了水文現(xiàn)象的物理成因,在原理上不同于分解合成法和還原或還現(xiàn)法。溯源重構(gòu)法采用源函數(shù)描述影響因子Xi作用于Y的規(guī)律。源函數(shù)是一種機理函數(shù),具有時不變的特性。同時溯源重構(gòu)法的基礎(chǔ)模型采用乘積模型,更符合水文系統(tǒng)的非線性相依特征。該方法物理含義明確,分析的不確定性相對較小。

    依據(jù)佳蘆河流域申家灣水文站的實測資料,應(yīng)用溯源重構(gòu)法實現(xiàn)了非一致性年最大洪峰流量序列向一致性的溯源重構(gòu)序列的轉(zhuǎn)化,結(jié)果表明溯源重構(gòu)法能夠有效消除年最大洪峰流量序列存在的突變及一階矩的線性或非線性趨勢,也能夠消除二階矩的非一致性。與分解合成法相比,該方法的目標也是試圖獲得一致性序列,但分解合成法僅能去除序列一階矩的線性趨勢,無法去除序列非線性趨勢以及方差的非一致性。相比于時變矩法確定出的時變概率分布,以及直接應(yīng)用非一致性原序列強行估計出的總體分布,溯源重構(gòu)序列的分布函數(shù)更符合原序列極值抽樣的統(tǒng)計屬性。因此,溯源重構(gòu)法所構(gòu)成的序列類似于原序列的模數(shù)序列,除可以用于水文頻率分析外,還可以通過對溯源重構(gòu)序列構(gòu)建預(yù)測模型,實現(xiàn)溯源重構(gòu)法在水文預(yù)測預(yù)報領(lǐng)域的應(yīng)用,并獲得不確定性相對較小的預(yù)報成果。事實上,只要科學合理的表達影響因子X,例如本文第3節(jié)中,若佳蘆河流域內(nèi)淤地壩發(fā)生漫壩或潰壩情形時,只需要對X,即有效產(chǎn)流面積,進行合理修正,就可以繼續(xù)沿用重構(gòu)函數(shù)獲得性質(zhì)相同的序列,因此溯源重構(gòu)法具有實際應(yīng)用前景。

    猜你喜歡
    洪峰流量水文一致性
    2022年《中國水文年報》發(fā)布
    關(guān)注減污降碳協(xié)同的一致性和整體性
    公民與法治(2022年5期)2022-07-29 00:47:28
    注重教、學、評一致性 提高一輪復(fù)習效率
    IOl-master 700和Pentacam測量Kappa角一致性分析
    水文
    水文水資源管理
    退耕還林工程對渭河洪峰流量的影響
    水文
    佛岡縣潖江流域年洪峰流量P-Ⅲ分布參數(shù)估算
    大南川流域設(shè)計洪峰流量計算分析
    日韩成人av中文字幕在线观看| 少妇人妻 视频| 久久女婷五月综合色啪小说| 亚洲人成77777在线视频| 免费少妇av软件| √禁漫天堂资源中文www| 国产av码专区亚洲av| 国产毛片在线视频| 久久国产精品男人的天堂亚洲 | 亚洲精品亚洲一区二区| 伦精品一区二区三区| 乱码一卡2卡4卡精品| 人妻系列 视频| 日韩强制内射视频| 国产精品一国产av| 欧美日韩成人在线一区二区| 国产女主播在线喷水免费视频网站| 亚洲不卡免费看| 日韩 亚洲 欧美在线| 免费观看的影片在线观看| 免费久久久久久久精品成人欧美视频 | 最新中文字幕久久久久| 国产一区二区在线观看日韩| 国产成人精品无人区| 婷婷成人精品国产| 男女高潮啪啪啪动态图| 26uuu在线亚洲综合色| 亚洲第一区二区三区不卡| 精品熟女少妇av免费看| 欧美最新免费一区二区三区| 中文字幕人妻熟人妻熟丝袜美| 麻豆精品久久久久久蜜桃| 中文字幕av电影在线播放| 亚洲天堂av无毛| 18在线观看网站| av视频免费观看在线观看| 久久人人爽人人片av| a 毛片基地| 国产精品一区www在线观看| 中文字幕最新亚洲高清| 国产熟女午夜一区二区三区 | 国产 一区精品| 日日爽夜夜爽网站| 母亲3免费完整高清在线观看 | 国产精品一国产av| 亚洲欧美一区二区三区黑人 | 国产高清国产精品国产三级| 美女国产视频在线观看| 男女边摸边吃奶| 精品国产一区二区三区久久久樱花| 亚洲精华国产精华液的使用体验| 久久精品国产亚洲av天美| 九色亚洲精品在线播放| 热re99久久国产66热| 日韩强制内射视频| 99热国产这里只有精品6| 日韩av不卡免费在线播放| 欧美最新免费一区二区三区| 国产精品一区二区在线观看99| 亚洲经典国产精华液单| 青春草视频在线免费观看| 久久精品夜色国产| 国产成人精品婷婷| 国产一区二区在线观看av| 日韩中文字幕视频在线看片| 黑人欧美特级aaaaaa片| 午夜老司机福利剧场| 国产精品一区二区在线不卡| 亚洲精品乱久久久久久| 高清av免费在线| 久久久久视频综合| 久久久久人妻精品一区果冻| 国产成人a∨麻豆精品| 一级毛片我不卡| 97精品久久久久久久久久精品| 少妇猛男粗大的猛烈进出视频| av播播在线观看一区| 久久免费观看电影| 满18在线观看网站| 国产精品无大码| 一区二区三区精品91| 国产在线免费精品| 日韩欧美一区视频在线观看| 午夜福利视频在线观看免费| 蜜桃久久精品国产亚洲av| 妹子高潮喷水视频| 秋霞在线观看毛片| 水蜜桃什么品种好| 国内精品宾馆在线| 99热全是精品| 十八禁高潮呻吟视频| 日韩精品免费视频一区二区三区 | 午夜激情久久久久久久| 成人国语在线视频| 亚洲国产欧美在线一区| 日韩精品有码人妻一区| 国产国语露脸激情在线看| 免费不卡的大黄色大毛片视频在线观看| 国产国拍精品亚洲av在线观看| 亚洲情色 制服丝袜| 18禁在线播放成人免费| 亚洲精品乱码久久久久久按摩| 国产乱来视频区| 日本欧美国产在线视频| 18禁动态无遮挡网站| 在线观看一区二区三区激情| 中国美白少妇内射xxxbb| 日日啪夜夜爽| videos熟女内射| 国产精品麻豆人妻色哟哟久久| av黄色大香蕉| 狂野欧美白嫩少妇大欣赏| 建设人人有责人人尽责人人享有的| 女性被躁到高潮视频| 91在线精品国自产拍蜜月| 美女xxoo啪啪120秒动态图| 美女福利国产在线| 国产亚洲一区二区精品| 精品99又大又爽又粗少妇毛片| 日本欧美国产在线视频| 精品少妇内射三级| 久久99热这里只频精品6学生| 91久久精品国产一区二区三区| 一边亲一边摸免费视频| 18禁裸乳无遮挡动漫免费视频| 搡女人真爽免费视频火全软件| 男人操女人黄网站| 久久99蜜桃精品久久| 亚洲精品视频女| 久久久亚洲精品成人影院| 国产毛片在线视频| 免费少妇av软件| 免费黄频网站在线观看国产| 午夜福利在线观看免费完整高清在| 国产免费又黄又爽又色| 伊人久久国产一区二区| 亚洲美女搞黄在线观看| 如何舔出高潮| 中文字幕制服av| 国产精品国产三级国产专区5o| 狂野欧美激情性bbbbbb| av女优亚洲男人天堂| 飞空精品影院首页| 精品视频人人做人人爽| 两个人免费观看高清视频| 欧美精品高潮呻吟av久久| 老女人水多毛片| 精品人妻在线不人妻| 久久女婷五月综合色啪小说| 天天操日日干夜夜撸| av免费观看日本| 五月玫瑰六月丁香| 日韩制服骚丝袜av| 3wmmmm亚洲av在线观看| 一边亲一边摸免费视频| 99九九在线精品视频| 国产精品蜜桃在线观看| 日韩一区二区视频免费看| 亚洲精品一二三| 亚洲欧美成人精品一区二区| 一级毛片aaaaaa免费看小| 久久毛片免费看一区二区三区| 欧美日韩一区二区视频在线观看视频在线| 亚洲av日韩在线播放| 国产极品天堂在线| 伊人亚洲综合成人网| a级毛片在线看网站| av在线观看视频网站免费| 超色免费av| 久久精品国产自在天天线| 日韩伦理黄色片| 极品人妻少妇av视频| 草草在线视频免费看| 久久久久久伊人网av| 天堂8中文在线网| 成年人午夜在线观看视频| 免费黄网站久久成人精品| 亚洲无线观看免费| 久久国产精品大桥未久av| 久久久久久久久久久免费av| 欧美激情极品国产一区二区三区 | 久久久久久久久久久丰满| 精品视频人人做人人爽| 全区人妻精品视频| 欧美bdsm另类| 丰满迷人的少妇在线观看| 亚洲国产色片| 亚洲欧美色中文字幕在线| 新久久久久国产一级毛片| 自线自在国产av| 午夜精品国产一区二区电影| 亚洲精品一区蜜桃| 日本欧美视频一区| 久久国产精品大桥未久av| 亚洲熟女精品中文字幕| 啦啦啦在线观看免费高清www| 久久人人爽人人爽人人片va| 国产视频首页在线观看| 亚洲av免费高清在线观看| 下体分泌物呈黄色| 性色av一级| 久久久久久久亚洲中文字幕| 久久狼人影院| 高清视频免费观看一区二区| 两个人的视频大全免费| 亚洲av福利一区| 久久午夜福利片| 乱码一卡2卡4卡精品| 少妇被粗大的猛进出69影院 | 久久 成人 亚洲| 成人无遮挡网站| 精品一区二区三区视频在线| 国产乱来视频区| 亚洲人成网站在线播| 国产日韩欧美亚洲二区| 精品一品国产午夜福利视频| 国产成人一区二区在线| 久久精品人人爽人人爽视色| 国产精品久久久久久久电影| 22中文网久久字幕| 午夜av观看不卡| 99久久中文字幕三级久久日本| 精品久久国产蜜桃| 日韩成人av中文字幕在线观看| 大片免费播放器 马上看| 五月天丁香电影| 性色av一级| 久久久久久人妻| 黄色欧美视频在线观看| 欧美激情国产日韩精品一区| 精品国产一区二区三区久久久樱花| 亚洲成人av在线免费| 男的添女的下面高潮视频| 日韩不卡一区二区三区视频在线| 久久久精品94久久精品| 日韩强制内射视频| 黄片播放在线免费| 日日爽夜夜爽网站| 老司机影院成人| 欧美精品一区二区免费开放| 亚洲av成人精品一二三区| 久久精品久久精品一区二区三区| 在线免费观看不下载黄p国产| 夜夜爽夜夜爽视频| 99re6热这里在线精品视频| 国产欧美亚洲国产| 一区二区三区精品91| 一级毛片电影观看| 99久国产av精品国产电影| 亚洲国产毛片av蜜桃av| 99国产精品免费福利视频| 大话2 男鬼变身卡| av网站免费在线观看视频| 中文欧美无线码| 国产男人的电影天堂91| 国产片特级美女逼逼视频| 午夜福利视频在线观看免费| 最近中文字幕2019免费版| 天堂中文最新版在线下载| 边亲边吃奶的免费视频| 69精品国产乱码久久久| 国产精品三级大全| 一级毛片我不卡| 亚洲国产欧美日韩在线播放| 午夜视频国产福利| 一级毛片黄色毛片免费观看视频| 这个男人来自地球电影免费观看 | 亚洲怡红院男人天堂| 男男h啪啪无遮挡| 老熟女久久久| 日韩视频在线欧美| 精品少妇黑人巨大在线播放| 亚洲欧美一区二区三区黑人 | 色婷婷久久久亚洲欧美| 欧美成人午夜免费资源| 日韩熟女老妇一区二区性免费视频| 精品国产乱码久久久久久小说| 最近中文字幕高清免费大全6| 一本大道久久a久久精品| 国模一区二区三区四区视频| 老司机亚洲免费影院| 妹子高潮喷水视频| 久久这里有精品视频免费| 午夜福利网站1000一区二区三区| 国产精品人妻久久久影院| 成年人午夜在线观看视频| 国产日韩欧美亚洲二区| 18+在线观看网站| 两个人的视频大全免费| 男女边吃奶边做爰视频| 性色avwww在线观看| 久久久久久久久久久丰满| 欧美激情极品国产一区二区三区 | 成人亚洲欧美一区二区av| 亚洲色图 男人天堂 中文字幕 | 超色免费av| 9色porny在线观看| 亚洲精品日本国产第一区| 99热全是精品| 欧美97在线视频| 亚洲成色77777| 欧美日韩在线观看h| 777米奇影视久久| 好男人视频免费观看在线| 久久久精品区二区三区| 大又大粗又爽又黄少妇毛片口| 人人妻人人澡人人爽人人夜夜| 日本与韩国留学比较| 国产国拍精品亚洲av在线观看| 一区二区三区免费毛片| 国产精品一二三区在线看| 男女无遮挡免费网站观看| 九九在线视频观看精品| 26uuu在线亚洲综合色| 国产日韩一区二区三区精品不卡 | 亚洲国产精品国产精品| 欧美激情国产日韩精品一区| 男女啪啪激烈高潮av片| 色吧在线观看| 国产成人免费无遮挡视频| 亚洲欧美日韩另类电影网站| 一个人免费看片子| 人妻夜夜爽99麻豆av| 精品久久蜜臀av无| 日韩 亚洲 欧美在线| av在线老鸭窝| 免费大片黄手机在线观看| av在线观看视频网站免费| 制服人妻中文乱码| 国产亚洲一区二区精品| 欧美+日韩+精品| 如日韩欧美国产精品一区二区三区 | 国产伦精品一区二区三区视频9| 免费观看无遮挡的男女| 国产在线一区二区三区精| 亚洲第一区二区三区不卡| 亚洲国产av新网站| 亚洲,一卡二卡三卡| 如何舔出高潮| 99久久精品一区二区三区| 一级,二级,三级黄色视频| 中国三级夫妇交换| xxxhd国产人妻xxx| 久热久热在线精品观看| 欧美日本中文国产一区发布| 免费观看a级毛片全部| 18禁在线无遮挡免费观看视频| 国产欧美日韩一区二区三区在线 | 十分钟在线观看高清视频www| 国产成人av激情在线播放 | 韩国高清视频一区二区三区| 亚洲怡红院男人天堂| 美女内射精品一级片tv| 国产深夜福利视频在线观看| 久久女婷五月综合色啪小说| 亚洲欧洲国产日韩| 久久国产亚洲av麻豆专区| 在线亚洲精品国产二区图片欧美 | 精品久久国产蜜桃| 欧美丝袜亚洲另类| 国产在线免费精品| 一区二区av电影网| av国产久精品久网站免费入址| 十八禁高潮呻吟视频| 午夜福利网站1000一区二区三区| 久久久精品94久久精品| 久久99热6这里只有精品| 日韩成人伦理影院| 日本与韩国留学比较| 亚洲精品av麻豆狂野| 免费黄频网站在线观看国产| 日韩av不卡免费在线播放| 国产综合精华液| 下体分泌物呈黄色| 久久精品熟女亚洲av麻豆精品| 成人亚洲欧美一区二区av| 亚洲欧洲国产日韩| 色婷婷久久久亚洲欧美| 国产精品99久久久久久久久| 女人精品久久久久毛片| 熟女人妻精品中文字幕| 欧美激情国产日韩精品一区| 一级毛片 在线播放| 欧美精品亚洲一区二区| 赤兔流量卡办理| 婷婷色综合大香蕉| 亚洲婷婷狠狠爱综合网| 18+在线观看网站| 少妇被粗大猛烈的视频| 亚洲精品av麻豆狂野| 欧美精品高潮呻吟av久久| 免费人妻精品一区二区三区视频| 午夜老司机福利剧场| 亚洲三级黄色毛片| 婷婷色麻豆天堂久久| 美女福利国产在线| 成人综合一区亚洲| 欧美日韩视频精品一区| 我要看黄色一级片免费的| 国产欧美亚洲国产| 久久久国产一区二区| 99热这里只有精品一区| 在线看a的网站| 国产日韩一区二区三区精品不卡 | 免费黄频网站在线观看国产| 满18在线观看网站| 黑人猛操日本美女一级片| 国产视频内射| 九草在线视频观看| 日日撸夜夜添| a级毛片黄视频| 欧美老熟妇乱子伦牲交| 哪个播放器可以免费观看大片| 亚洲精品av麻豆狂野| 国产亚洲精品第一综合不卡 | 国产av码专区亚洲av| 777米奇影视久久| 少妇人妻精品综合一区二区| 赤兔流量卡办理| 久久久久精品久久久久真实原创| 大片免费播放器 马上看| 99久久精品一区二区三区| 3wmmmm亚洲av在线观看| 熟女人妻精品中文字幕| 亚洲经典国产精华液单| 亚洲欧洲精品一区二区精品久久久 | 精品熟女少妇av免费看| 国产高清国产精品国产三级| 国产深夜福利视频在线观看| 视频在线观看一区二区三区| 丝袜美足系列| 国产精品不卡视频一区二区| 国产精品免费大片| 伊人久久国产一区二区| 一级,二级,三级黄色视频| 黄色毛片三级朝国网站| 成人综合一区亚洲| 欧美bdsm另类| 久久久久国产精品人妻一区二区| 成人亚洲欧美一区二区av| 校园人妻丝袜中文字幕| 久久久精品免费免费高清| 亚洲国产精品专区欧美| 亚洲欧美一区二区三区国产| 尾随美女入室| 亚洲av.av天堂| 最近2019中文字幕mv第一页| 欧美+日韩+精品| 人成视频在线观看免费观看| 综合色丁香网| 成年女人在线观看亚洲视频| 欧美亚洲 丝袜 人妻 在线| √禁漫天堂资源中文www| 香蕉精品网在线| www.av在线官网国产| 看非洲黑人一级黄片| 国产日韩欧美视频二区| 亚洲av欧美aⅴ国产| 午夜老司机福利剧场| 国产在线视频一区二区| 精品一区在线观看国产| 亚洲第一av免费看| 最黄视频免费看| 欧美精品亚洲一区二区| 国产成人精品在线电影| 大片电影免费在线观看免费| 久久久久久人妻| 一区在线观看完整版| 人妻系列 视频| 高清黄色对白视频在线免费看| 黄片播放在线免费| 黄色欧美视频在线观看| 成人午夜精彩视频在线观看| 26uuu在线亚洲综合色| 国产乱来视频区| 亚洲精品aⅴ在线观看| 能在线免费看毛片的网站| 人人妻人人爽人人添夜夜欢视频| 爱豆传媒免费全集在线观看| 中文字幕制服av| av在线观看视频网站免费| 久久久久网色| 老司机影院毛片| 亚洲av电影在线观看一区二区三区| 在线观看人妻少妇| 欧美bdsm另类| 韩国高清视频一区二区三区| 在线看a的网站| 卡戴珊不雅视频在线播放| 少妇的逼水好多| 中文精品一卡2卡3卡4更新| 丝瓜视频免费看黄片| 国产欧美日韩综合在线一区二区| 天堂俺去俺来也www色官网| 久久久午夜欧美精品| 老司机影院成人| 欧美三级亚洲精品| 热re99久久国产66热| a级片在线免费高清观看视频| 亚洲av.av天堂| 啦啦啦啦在线视频资源| 人妻夜夜爽99麻豆av| 91精品伊人久久大香线蕉| 两个人免费观看高清视频| 国产成人91sexporn| 免费看av在线观看网站| 国产成人精品无人区| 毛片一级片免费看久久久久| 国产精品人妻久久久久久| 成人国产av品久久久| 一级爰片在线观看| 久久精品国产自在天天线| 老司机影院成人| 国产白丝娇喘喷水9色精品| 制服人妻中文乱码| 国产成人精品一,二区| 最黄视频免费看| 精品人妻熟女av久视频| 国产高清不卡午夜福利| 成人漫画全彩无遮挡| 中文字幕久久专区| 99热国产这里只有精品6| 国产精品一区www在线观看| 久久人妻熟女aⅴ| 少妇猛男粗大的猛烈进出视频| 91精品三级在线观看| 国产高清国产精品国产三级| 秋霞伦理黄片| 又粗又硬又长又爽又黄的视频| 久久久a久久爽久久v久久| 午夜av观看不卡| 水蜜桃什么品种好| 肉色欧美久久久久久久蜜桃| 午夜福利视频精品| 看免费成人av毛片| 国产成人av激情在线播放 | .国产精品久久| 久久午夜福利片| 一级毛片黄色毛片免费观看视频| 久久 成人 亚洲| 欧美激情极品国产一区二区三区 | 久久久精品94久久精品| 在线观看免费视频网站a站| 2018国产大陆天天弄谢| 精品久久久噜噜| 精品午夜福利在线看| 亚洲av成人精品一区久久| 各种免费的搞黄视频| 久久国产精品男人的天堂亚洲 | 国产成人91sexporn| 九九爱精品视频在线观看| 91aial.com中文字幕在线观看| 51国产日韩欧美| 91久久精品电影网| 国产精品熟女久久久久浪| 欧美xxⅹ黑人| 免费日韩欧美在线观看| 日本-黄色视频高清免费观看| 能在线免费看毛片的网站| 免费高清在线观看日韩| 人人澡人人妻人| 午夜视频国产福利| 欧美精品国产亚洲| 老女人水多毛片| 七月丁香在线播放| 乱人伦中国视频| 亚洲欧洲国产日韩| 亚洲不卡免费看| 国产精品秋霞免费鲁丝片| 国产欧美日韩一区二区三区在线 | 曰老女人黄片| 三级国产精品欧美在线观看| 精品卡一卡二卡四卡免费| 黄色配什么色好看| 观看美女的网站| 欧美日本中文国产一区发布| 高清不卡的av网站| 免费av不卡在线播放| 最后的刺客免费高清国语| 全区人妻精品视频| 亚洲精品日本国产第一区| 亚洲av欧美aⅴ国产| a 毛片基地| 亚洲精品视频女| 久热久热在线精品观看| 少妇人妻久久综合中文| 精品人妻一区二区三区麻豆| 免费播放大片免费观看视频在线观看| 美女视频免费永久观看网站| 26uuu在线亚洲综合色| 乱人伦中国视频| 97超视频在线观看视频| 久久久久久久亚洲中文字幕| 精品午夜福利在线看| 亚洲欧美色中文字幕在线| 日韩中文字幕视频在线看片| 美女中出高潮动态图| av天堂久久9| 亚洲久久久国产精品| 啦啦啦视频在线资源免费观看| 色94色欧美一区二区| 夜夜骑夜夜射夜夜干| 亚洲精品aⅴ在线观看| 亚洲欧美一区二区三区国产| 久久免费观看电影| 一级黄片播放器| 汤姆久久久久久久影院中文字幕| a 毛片基地| 国产黄色免费在线视频| 日韩一区二区三区影片| 日本与韩国留学比较| 国产一级毛片在线| 搡女人真爽免费视频火全软件| 又黄又爽又刺激的免费视频.| 精品久久久久久久久av|