秦 毅,李 時
(省部共建西北旱區(qū)生態(tài)水利國家重點實驗室西安理工大學(xué),陜西西安 710048)
近年來,大量研究成果表明實測水文序列存在明顯的非一致性特征[1-4]。它帶來兩方面問題:(1)在非一致性條件下,繼續(xù)使用以一致性樣本為基礎(chǔ)的預(yù)測分析方法,包括傳統(tǒng)水文頻率分析方法、多元統(tǒng)計分析的Coupla方法、機(jī)器學(xué)習(xí)等方法是否能獲得合理的成果,以及對未來的預(yù)測應(yīng)該如何評價。(2)水文模型參數(shù)的率定工作變得困難,因為參數(shù)率定采用資料的非一致性導(dǎo)致了不確定性的增加。相比較而言,第一個問題在設(shè)計分析方面顯得更為突出,例如水利工程、水土保持工程及生態(tài)工程的設(shè)計,洪旱災(zāi)害的防治等都離不開對實測樣本序列的統(tǒng)計分析與預(yù)測,除非研究對象的物理機(jī)制被真實地掌握,否則所做出的水文分析與預(yù)測的結(jié)果只有在一致性(平穩(wěn)性)序列條件下才有意義。因此學(xué)者們投入了大量研究,包括:(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)方法外,近年來有學(xué)者們將研究焦點放在了時變矩法[11-12]上,并據(jù)此建立時變概率分布,再通過重現(xiàn)期或可靠度與時變概率分布間的關(guān)系確定出工程設(shè)計值。當(dāng)前比較有代表性的方法包括基于修正重現(xiàn)期的EWT[13-14]和ENE方法[15-16],和基于可靠度的ER[17]、DLL[18]和ADLL[19]等方法。
然而,這些方法仍然存在一定的爭議。首先是對趨勢的看法,有學(xué)者認(rèn)為觀測到的水文序列的趨勢并非是非一致性特征的表現(xiàn),而可能是平穩(wěn)過程中的周期性頻率波動所致[20],甚至認(rèn)為趨勢這個詞都缺少明確的定義[21];再者,實踐中,在工程設(shè)計壽命內(nèi),一定可靠度下的設(shè)計值將隨著計算時段的起始時間以及統(tǒng)計參數(shù)與協(xié)變量間關(guān)系擬合曲線的類型而變化,這就意味著未來設(shè)計值的可靠性嚴(yán)重依賴統(tǒng)計參數(shù)的時變特征,而這種時變特征卻因為時間序列遍歷性的缺失,造成對未來水文序列統(tǒng)計參數(shù)預(yù)測的不確定性[22]。正如Serinaldi等[23]指出的,如果模型結(jié)構(gòu)無法通過演繹推斷的方式獲得,特別是描述分布參數(shù)非平穩(wěn)變化模型是僅通過歸納推理擬合得來,則一個額外不確定性來源被無形中引入了模型結(jié)構(gòu)中,使得所得到的非平穩(wěn)模型無法在實際應(yīng)用中提高設(shè)計值的可信度與準(zhǔn)確性,而模型的這種不確定性將很容易導(dǎo)致不具物理意義的結(jié)果產(chǎn)生;此外,利用時變參數(shù)識別非平穩(wěn)序列的分布函數(shù)過程中,所采用的經(jīng)驗頻率計算方法也是引起不確定性的重要因素,因為計算經(jīng)驗頻率的方法也同樣要求樣本為簡單隨機(jī)樣本,基于非平穩(wěn)序列直接進(jìn)行經(jīng)驗頻率計算顯然不恰當(dāng)??梢姡瑢τ诮y(tǒng)計特征隨時間變化的單一隨機(jī)變量,其分布研究尚存在如此多的問題,更何況解決多元時變隨機(jī)變量的聯(lián)合分布問題,難度可想而知。
很顯然,造成“難度”的核心是樣本序列的非簡單性。如果能夠?qū)⒎且恢滦詷颖拘蛄修D(zhuǎn)換為一致性序列,則因為對此簡單序列已經(jīng)有了成熟的分析理論、技術(shù)和用好技術(shù)的方法,前述困難將不復(fù)存在。因此,本文提出溯源重構(gòu)法,它是基于反映水文現(xiàn)象物理機(jī)理的源函數(shù)將非一致性序列重構(gòu)為新的一致性序列的方法。依靠重構(gòu)的一致性序列,解決需要基于一致性序列才能進(jìn)行分析工作的問題。事實上,當(dāng)前已存在將非一致性序列轉(zhuǎn)變?yōu)橐恢滦孕蛄械姆椒?,即還原或還現(xiàn)法、分解合成法,本文將在后面就這兩種方法與溯源重構(gòu)法進(jìn)行比較與討論。
2.1 源函數(shù)的概念一般來說,一個水文變量Y的變化是多個影響因子X1,…,XN共同作用的結(jié)果,見圖1。在排除其他影響因子影響的條件下,單獨考察影響因子Xi對水文變量Y的作用,將這種作用規(guī)律記為fi(Xi),則fi(Xi)是Xi對Y的物理作用機(jī)制,因此這里稱fi(Xi)是Xi作用于Y的源函數(shù)。作為物理機(jī)制,源函數(shù)是不會改變的,就像河道流量Q=AV,其中A是過流斷面的面積,V是流速,則V和A各自對Q的作用規(guī)律或說機(jī)制,也即源函數(shù)就是fV(V)=V,fA(A)=A。當(dāng)V或A或兩者共同隨時間變化時,流量也會隨著V(t)或A(t)或兩者在不同時刻的取值而表現(xiàn)出隨時間變化的Q(t)。但是,不論V和A怎樣取值,源函數(shù)fV(V)=V,fA(A)=A保持不變,而流量等于這兩個源函數(shù)相乘的機(jī)制也不會改變。即物理機(jī)制不隨時間變化,變化的永遠(yuǎn)是Xi的狀態(tài)。
圖1 各影響因子對研究變量Y的源函數(shù)作用關(guān)系
嚴(yán)格的說,源函數(shù)的確定有兩種方法:理論推導(dǎo)(數(shù)學(xué)解析)和實驗分析。實際上,由于水文現(xiàn)象變化的復(fù)雜性,即便我們認(rèn)可源函數(shù)的存在,也會因為人類對自然的有限認(rèn)知而很難得到絕對準(zhǔn)確的物理機(jī)制數(shù)學(xué)表達(dá),至少當(dāng)前很少看到對單一因子影響的研究。由于事物之間的作用規(guī)律可以隱含于統(tǒng)計規(guī)律中,所以當(dāng)下可以借用Y和Xi間的統(tǒng)計關(guān)系來近似估計源函數(shù)。
既然源函數(shù)或由源函數(shù)構(gòu)成的函數(shù)不隨時間變化,我們就可以利用這個特征將非一致性序列轉(zhuǎn)變?yōu)橐恢滦孕蛄小?/p>
2.2 溯源重構(gòu)方法的基礎(chǔ)模型欲利用源函數(shù)時不變的特點,就要表達(dá)出多個影響因子各自對Y作用后的一般綜合關(guān)系。鑒于數(shù)學(xué)解析方法暫時行不通,故采用建立統(tǒng)計模型的方法來表征研究對象Y與解釋變量源函數(shù)之間的關(guān)系。眾所周知,數(shù)學(xué)的基本運算包括加(減)與乘(除),結(jié)合運算法則,Y與解釋變量源函數(shù)之間的關(guān)系有兩種一般性廣義表達(dá)方式,一是各解釋變量間相互獨立時采用的疊加模式:
另一個是可以考慮解釋變量間存在相關(guān)性的乘積模式:
其中:N指作用于研究變量Y的所有影響因子的個數(shù),這些影響因子包括人類認(rèn)識到的和尚未認(rèn)識的;fi(Xi(t))為影響因子Xi的源函數(shù)。假設(shè)在N個影響因子中,有n個因子是人類目前已認(rè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(rèn)為是自然隨機(jī)變量,具有概率分布,并具有對應(yīng)的均值μ、μ*,和方差σ2、σ*2。令Ks(t)=f1(X1(t))+…+fn(Xn(t))和Kp(t)=f1(X1(t))…fn(Xn(t))。在生產(chǎn)實踐當(dāng)中,由于我們更關(guān)心影響因子取給定值時,Y的變化情況,與回歸分析的考慮相同,在計算Y之前,解釋變量X1(t),…,Xn(t)的取值已經(jīng)確定,可以作為可控變量處理[24]。于是Ks(t)和Kp(t)在t時刻為常數(shù),故Ks(t)和Kp(t)是已知的時間函數(shù)。重新寫上述式(3)、式(4)為:
對比兩種表達(dá)式下Y(t)的條件數(shù)學(xué)期望和條件方差如下。
疊加模型:
乘積模型:
不難看到,兩種模型下Y(t)的條件均值均是時變函數(shù),一旦Xi發(fā)生趨勢性變化,Y(t)的條件數(shù)學(xué)期望也就出現(xiàn)趨勢性變化;與均值不同,在疊加模型下的方差為常數(shù),而乘積模型下方差是時變函數(shù),且會隨Xi的趨勢性變化而產(chǎn)生更強(qiáng)的趨勢性。
現(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],特別是極端事件的極端特性越強(qiáng),二階矩的變化就越不能被忽視[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)剔除全部非一致性因子影響后剩余一致性影響因子以及尚未被人類認(rèn)知的影響因子源函數(shù)的乘積,它是平穩(wěn)的隨機(jī)過程。如果只剔除了部分非一致性變化因子X1,X2,…,Xl(l 理論上,當(dāng)引起Y非一致性變化的影響因子全部被識別出來,即l=m,且這些因子的源函數(shù)fi(Xi(t))已有準(zhǔn)確的表達(dá),那么由溯源重構(gòu)函數(shù)計算出來的溯源重構(gòu)序列RSt={rs1,rs2,…,rsn}為一致性序列,可用于需要一致性序列才能進(jìn)行工作的情況。生產(chǎn)實踐中,因人類認(rèn)知和現(xiàn)有方法的局限,我們不可能得到絕對的平穩(wěn)過程,這也不是我們的目標(biāo),但隨著非一致性變化的解釋變量對應(yīng)的正確源函數(shù)不斷加入溯源重構(gòu)函數(shù)中,則溯源重構(gòu)序列將逐漸趨向并最終趨向于平穩(wěn),接近于一個隨機(jī)噪聲。 需要提醒的是,由于估計源函數(shù)存在不確定性,所以使用中必須配合溯源重構(gòu)序列RSt的一致性檢驗。若檢驗結(jié)果為一致,則同時接受所確定的影響因子及其源函數(shù)的估計形式,否則,仍需進(jìn)一步研究影響因子及其源函數(shù)。 由于研究變量Y和影響因子間的因果關(guān)系,研究變量Y的各種非一致性變化(線性、非線性、突變)與影響因子的相應(yīng)非一致性變化一致,這使得溯源重構(gòu)法去除非一致性時并不受非一致性變化類型的影響,即具有普遍性。相較于對非一致性變化的Y進(jìn)行純數(shù)學(xué)模擬的時變矩法,溯源重構(gòu)方法更注重找到引起Y變化的非一致性變化影響因子Xi和它對Y的作用規(guī)律,并以此規(guī)律將非一致性變化從Y序列中排除,這正是溯源重構(gòu)法能夠?qū)⒎且恢滦宰兓腨(t)變?yōu)橐恢滦孕蛄蠷S(t)的關(guān)鍵所在。 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)研等手段,初步確定研究變量的主要影響因子,之后對這些主要影響因子進(jìn)行趨勢分析,找出造成水文變量非一致性變化的影響因子。 佳蘆河位于黃河中游的河龍區(qū)間,毛烏素沙漠南緣,陜北黃土高原北段,屬于黃土高原溝壑區(qū)。河流由西北流向東南,至佳縣佳蘆鎮(zhèn)木廠灣村注入黃河。河流全長93 km,流域面積1134 km2,出口控制站為申家灣水文站。流域洪水由降雨引發(fā)。由于氣候較為干旱,自然條件較差,該地區(qū)容易發(fā)生風(fēng)沙災(zāi)害,水土流失現(xiàn)象普遍。為了控制水土流失,流域內(nèi)修建了數(shù)十座淤地壩。淤地壩一般由兩大件組成——壩和泄洪臥管。為了使入庫泥沙得到沉積,實現(xiàn)多攔沙少攔水的功效,淤地壩運行方式設(shè)計成:發(fā)生洪水時,高含沙洪水首先被攔截在庫內(nèi),后由臥管分?jǐn)?shù)天將洪水排出庫外,且要求臥管排水速度緩慢,排水流量一般為1~1.5 m3/s。結(jié)果導(dǎo)致出口斷面洪峰較建壩前明顯減小,庫壩泄流對流域出口的洪峰流量的貢獻(xiàn)可以忽略不計。這種效果等同于流域產(chǎn)生洪峰的產(chǎn)流面積減少,減少的數(shù)量為庫壩控制面積。這里稱流域面積與庫壩控制面積的差為有效產(chǎn)流面積;另一方面退耕還林以及淤地壩的蓄水,使得坡面較大范圍內(nèi)的土壤含水率增加,促進(jìn)了植被增長,導(dǎo)致流域植被截留、下滲、匯流等條件均發(fā)生改變,也造成了洪峰流量序列的非一致性變化。因此,選擇有效產(chǎn)流面積與植被覆蓋度作為佳蘆河流域年最大洪峰流量序列非一致性變化的主要影響因子。 洪峰流量序列采用申家灣水文站1957—2010年實測數(shù)據(jù)進(jìn)行分析。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)流面積的變化??梢钥吹?,流域洪峰流量和有效產(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ù)進(jìn)行分析。由于淤地壩因逐年修建而增多,且建壩本身不屬于隨機(jī)事件,所以有效產(chǎn)流面積序列呈現(xiàn)顯著單調(diào)下降的趨勢,也沒有一般隨機(jī)序列的波動變化。 圖3 佳蘆河流域有效產(chǎn)流面積序列 圖4 佳蘆河流域植被覆蓋度序列 檢驗結(jié)果(表1)表明年最大洪峰流量序列的均值和方差都出現(xiàn)了顯著減少趨勢,但植被覆蓋度的均值出現(xiàn)顯著上升趨勢。由于年最大洪峰流量一般在汛期發(fā)生,本文對佳蘆河流域申家灣站的各時段致洪降雨序列進(jìn)行了趨勢性檢驗??紤]到黃土丘陵溝壑區(qū)降雨過程一般歷時較短,不超過12 h,因此最大12 h降雨能夠代表相應(yīng)洪水的致洪暴雨,故在此處僅展示最大12 h降雨變化情況(圖2)。同時將與佳蘆河流域較近的榆林氣象站最大日降雨序列作為參考,一并進(jìn)行檢驗。檢驗結(jié)果(表2)顯示,佳蘆河流域致洪降雨序列未發(fā)生顯著變化(包括線性和非線性一階矩、二階矩)。因此排除了降雨變化造成洪峰流量非一致性變化的可能。 圖2 申家灣站年最大洪峰流量序列與致洪暴雨序列 表1 佳蘆河流域洪峰流量序列及其影響因子的趨勢性分析 表2 佳蘆河流域致洪降雨序列趨勢性分析 3.2 溯源重構(gòu)序列的生成通過對流域?qū)崪y資料的分析,我們可以得到這樣的認(rèn)識:佳蘆河流域的降雨未發(fā)生顯著的非一致性變化,申家灣站年最大洪峰流量序列的非一致性變化主要產(chǎn)生于流域有效產(chǎn)流面積和植被覆蓋度的顯著變化。因此,根據(jù)溯源重構(gòu)法的思想,重構(gòu)時不考慮降雨因子,僅考慮有效產(chǎn)流面積與植被覆蓋度兩個因子。在對洪峰流量序列進(jìn)行重構(gòu)之前,首先要分析出有效產(chǎn)流面積和植被覆蓋度對洪峰流量作用的源函數(shù)。 3.2.1 源函數(shù)的建立方法 源函數(shù)是溯源重構(gòu)法的關(guān)鍵,也是提升對水文過程認(rèn)知的重要基礎(chǔ),是不應(yīng)該被忽略的學(xué)術(shù)問題。它的確定應(yīng)該來自于理論推導(dǎo)或?qū)嶒?,但現(xiàn)階段,介于水文過程的復(fù)雜性,基于理論和實驗手段成功確定出的源函數(shù)幾乎看不到。因此,本著影響因子的作用規(guī)律必然體現(xiàn)在統(tǒng)計規(guī)律中的認(rèn)識,我們可以通過回歸的方法來逐個估計源函數(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)系估計,依次進(jìn)行,直到m個因子的源函數(shù)被估計出。 顯然,這種做法會帶來因影響因子的引入順序改變各因子源函數(shù)形式的問題。理論上,真實的源函數(shù)將不會受重構(gòu)順序的干擾。而實際中,由于源函數(shù)是估計函數(shù),重構(gòu)順序的干擾是可能存在的,為避免這個問題,可根據(jù)各影響因子與研究變量之間相關(guān)性的強(qiáng)弱順序確定重構(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)序列 當(dāng)確立了有效產(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)性的強(qiáng)弱順序確定重構(gòu)順序。將洪峰流量序列Qm,t及其影響因子序列Ae,t和FVCt的值分別代入式(18)—(20)中,得到相應(yīng)的溯源重構(gòu)序列RSA,t、RSFVC,t和RSdouble,t。 按照本文提及的溯源重構(gòu)理論建立的溯源重構(gòu)方法,其實質(zhì)是一種轉(zhuǎn)換方法,所獲得的重構(gòu)序列RSt理應(yīng)是一個具有一致性的隨機(jī)序列。作為驗證,本文以第3節(jié)中的溯源重構(gòu)序列RSA,t、RSFVC,t和RSdouble,t為對象,通過考察其趨勢性變化,探究溯源重構(gòu)理論的正確性和溯源重構(gòu)序列的合理性。同時,鑒于以下兩個原因,對重構(gòu)序列的特征進(jìn)行探討同樣顯得必要。①雖然RSt(如RSA,t、RSFVC,t和RSdouble,t)不是直接觀測變量,但因它是由觀測變量計算出來,因此有物理意義:相當(dāng)于模數(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)序列的分布特征分布特征是隨機(jī)變量取值特性的完整描述,是研究隨機(jī)現(xiàn)象的核心。這里討論年最大洪峰流量溯源重構(gòu)序列的分布函數(shù),并與基于時變矩獲得的洪峰流量序列(原序列)的分布函數(shù)進(jìn)行比較。同時,針對當(dāng)前工程師們無奈之選,即無視非平穩(wěn)的存在,直接進(jì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)度準(zhǔn)則確定最優(yōu)分布。申家灣站溯源重構(gòu)序列的最優(yōu)分布為GEV,對應(yīng)的分布參數(shù)分別為μ=9.657,σ=83.525,κ=0.878。對于非一致性序列的時變分布則基于GAMLSS模型[38]采用時變矩法確定,備選分布仍選取上述五種類型,以有效產(chǎn)流面積與植被覆蓋度作為協(xié)變量,評價準(zhǔn)則為AIC準(zhǔn)則[39],具體見表4。而無視非一致性存在,直接以非一致性原洪峰流量序列估計的總體分布也列于表4。 表4 不同方法確定的年最大洪峰流量序列的分布類型 理論上,年最大洪峰流量序列屬于極值序列,它的分布應(yīng)該是極值型。從表4可以看到不同屬性序列和不同分析方法所挑選出的分布函數(shù)不同,其中只有溯源重構(gòu)序列的分布符合極值型,本文作者對陜西省榆林地區(qū)大理河等6個流域的年最大洪峰流量溯源重構(gòu)序列進(jìn)行分布函數(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è)計值是水資源利用、管理與保護(hù)領(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)橐恢滦孕蛄校瑸橐砸恢滦孕蛄凶鳛榛A(chǔ)的分析計算工作提供樣本資料。這種想法同樣體現(xiàn)在還原或還現(xiàn)法和分解合成法中。但由于方法所依據(jù)的原理不同,方法的應(yīng)用效果隨之不同。 還原或還現(xiàn)法依據(jù)的原理是通過對水文變量系統(tǒng)性變化的增量回加到原序列或?qū)λ淖兞窟M(jìn)行修正,達(dá)到使水文變量序列一致的目的。在增量的量化或修正值的量化工作中,由于要求資料的種類多、數(shù)量大,且對資料的準(zhǔn)確性依賴程度高,使用不易。同時牽扯的一些概念難以精確界定[40],存在方法失真或失效[41]等問題。即便用分布式水文模型直接計算增量,也會面臨資料的非一致性導(dǎo)致模型參數(shù)率定(優(yōu)化法)的不確定性加大問題[41]。正因如此,本文不做具體的分析對比。 分解合成法[42]認(rèn)為,“無論水文現(xiàn)象的變化多么復(fù)雜,水文序列總可以分解成確定性成分和隨機(jī)性成分這兩部分。一般來說,水文序列的確定性成分主要受人類活動的影響,但并不排除氣候因素(如氣候轉(zhuǎn)型期)和下墊面因素(如火山爆發(fā)、地震等)的影響,其變化規(guī)律可以在較短的工程年代里發(fā)生緩慢的漸變或劇烈的突變,因此水文序列中確定性成分的變化規(guī)律往往是非一致的。而水文序列的隨機(jī)性成分則主要受氣候、地質(zhì)等因素的影響,其變化規(guī)律需要一個漫長的地質(zhì)年代才能改變,因此水文序列中隨機(jī)性成分的統(tǒng)計規(guī)律是相對一致的?!庇谑侵x平等[42]建議采用成因分析法與統(tǒng)計分析法分別對確定性成分和隨機(jī)性成分進(jìn)行識別與檢驗,并對確定性成分進(jìn)行擬合計算,對隨機(jī)性成分采用傳統(tǒng)方法進(jìn)行頻率分析計算。后者意味著分解合成法獲得的隨機(jī)變量樣本是一個簡單隨機(jī)樣本。下面將繼續(xù)以佳蘆河流域年最大洪峰流量序列資料,研究分析分解合成法生成的隨機(jī)樣本特征。 根據(jù)分解合成法原理,非一致性水文序列Xt可按下式表示: 式中:Yt為確定的非周期成分(包括趨勢、跳躍等);Pt為確定性的周期成分(包括簡單的或復(fù)合周期的成分等);St為隨機(jī)成分。 為避免周期對趨勢的影響,這里首先進(jìn)行周期分析,利用小波分析[43]得到周期序列存在21年、9年兩個顯著周期,由此獲得周期函數(shù)如式(22): 剔除周期后發(fā)現(xiàn)序列仍存在明顯的趨勢,通過回歸分析得到式(23)所示的趨勢成分: 根據(jù)式(21)剔除周期成分與趨勢成分后得到隨機(jī)序列St,見圖6,對該隨機(jī)序列進(jìn)行趨勢檢驗,結(jié)果見表5。 圖6 分解合成法構(gòu)建的佳蘆河年最大洪峰流量一致性序列 表5 分解合成法構(gòu)建的佳蘆河年最大洪峰流量序列的一致性檢驗 分解合成法獲得的佳蘆河年最大洪峰流量的重構(gòu)序列十分有效地消除了一階矩的顯著線性趨勢,但Spearman秩相關(guān)系數(shù)依然大于臨界值,表明其一階矩的非線性型非一致性仍然顯著存在;方差的非一致性檢驗統(tǒng)計量χ為23.064,幾乎沒有任何改善,說明分解合成法重構(gòu)的隨機(jī)成分并非簡單樣本。 產(chǎn)生上述現(xiàn)象的原因可以概化為兩個方面。一是由于定義兩種成分的物理機(jī)制均不夠清晰,例如分解合成法原理述及的“水文序列的確定性成分主要受人類活動的影響,但并不排除氣候因素(如氣候轉(zhuǎn)型期)和下墊面因素(如火山爆發(fā)、地震等)的影響”與“隨機(jī)性成分則主要受氣候、地質(zhì)等因素的影響,其變化規(guī)律需要一個漫長的地質(zhì)年代才能改變”,其中氣候轉(zhuǎn)型期以及需要漫長的地質(zhì)年代才能改變的氣候因素?zé)o法清晰界定,更不能明確表達(dá),因此確定性成分的擬合具有較大的不確定性,由實測資料推求出的隨機(jī)性成分為實測資料與確定性成分之差,就會具有較大的不確定性;二是該方法假定水文序列Xt的各個成分滿足線性疊加特性(即疊加模型)[44]。因此,根據(jù)上文疊加模型與乘積模型的對比分析,分解合成法構(gòu)建的“相對一致”性序列僅能實現(xiàn)均值一致,而方差非一致的現(xiàn)象就不難理解了。這也成為了疊加模型無法描述序列的二階矩變化的例證。 相比之下,溯源重構(gòu)法認(rèn)為水文變量Y的變化是一系列影響因子自身的非一致性變化引起的。這些因子的非一致性影響通過它們各自的源函數(shù)作用于Y。如果將全部非一致性變化因子的作用從Y序列中剔除,則剔除后的序列RSt就是一致性序列。與分解合成法不同,溯源重構(gòu)法中的隨機(jī)成分是指未被人類認(rèn)知的次要因子的作用,無需專門率定。同時溯源重構(gòu)法的基礎(chǔ)模型采用乘法模式,更符合水文系統(tǒng)的非線性相依特征。溯源重構(gòu)法的物理含義明確,分析的不確定性相對較小。 眾所周知,任何現(xiàn)象都有其因果關(guān)系,以數(shù)學(xué)術(shù)語來說,即因變量與自變量的關(guān)系。作者認(rèn)為,水文變量Y的非一致性變化是由于自變量(影響因子)自身的非一致性變化所引起的。據(jù)此提出溯源重構(gòu)法,即從源函數(shù)概念出發(fā),依據(jù)因變量與自變量間的物理機(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ù)是一種機(jī)理函數(shù),具有時不變的特性。同時溯源重構(gòu)法的基礎(chǔ)模型采用乘積模型,更符合水文系統(tǒng)的非線性相依特征。該方法物理含義明確,分析的不確定性相對較小。 依據(jù)佳蘆河流域申家灣水文站的實測資料,應(yīng)用溯源重構(gòu)法實現(xiàn)了非一致性年最大洪峰流量序列向一致性的溯源重構(gòu)序列的轉(zhuǎn)化,結(jié)果表明溯源重構(gòu)法能夠有效消除年最大洪峰流量序列存在的突變及一階矩的線性或非線性趨勢,也能夠消除二階矩的非一致性。與分解合成法相比,該方法的目標(biāo)也是試圖獲得一致性序列,但分解合成法僅能去除序列一階矩的線性趨勢,無法去除序列非線性趨勢以及方差的非一致性。相比于時變矩法確定出的時變概率分布,以及直接應(yīng)用非一致性原序列強(qiá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ù)報成果。事實上,只要科學(xué)合理的表達(dá)影響因子X,例如本文第3節(jié)中,若佳蘆河流域內(nèi)淤地壩發(fā)生漫壩或潰壩情形時,只需要對X,即有效產(chǎn)流面積,進(jìn)行合理修正,就可以繼續(xù)沿用重構(gòu)函數(shù)獲得性質(zhì)相同的序列,因此溯源重構(gòu)法具有實際應(yīng)用前景。3 實例研究
4 討論
5 結(jié)論