本發(fā)明涉及一種月徑流模擬方法，具體的說是針對全河任意點的月徑流時間位移二維矩陣的模擬方法。

背景技術：

流域水資源統(tǒng)一調(diào)度，首先需要厘清復雜水文系統(tǒng)的風險及不確定性，水文系統(tǒng)的不確定性，包括隨機性、模糊性、灰色性和未確定性，一直以來都是水文研究的熱點。河流來水過程對水庫調(diào)度起著至關重要的影響，其隨機性正是水文不確定的一種。隨機過程的理論和時間序列分析技術可以用來研究徑流過程的統(tǒng)計變化特性，為水文水資源系統(tǒng)的規(guī)劃、調(diào)度提供理論依據(jù)，是認識、設計和管理復雜水文系統(tǒng)的主要方式之一。

徑流模擬序列作為水資源統(tǒng)一調(diào)度的輸入條件，不僅為水庫群提供來水信息，還需為沿河各取水口點提供供水信息，而這些點中只有少部分有徑流數(shù)據(jù)資料。現(xiàn)有的多站點徑流模擬模型均是依靠現(xiàn)有歷史資料對已建水文站點進行徑流模擬，已不能滿足水資源統(tǒng)一調(diào)度的需求，這就需要開發(fā)一種新的徑流模擬模型，通過分析已建站點的自相關關系和站點間互相關關系，采用一定的方法進行降階，從而模擬全河任意點的徑流時間序列。

技術實現(xiàn)要素：

本發(fā)明所要解決的技術問題是克服現(xiàn)有技術的缺陷，提供一種適用于當前氣候變化條件下非線性、偏態(tài)分布的時間位移二維月徑流模擬模型，該模型不僅可模擬單站點的月徑流時間序列、多站點的月徑流聯(lián)合分布，還可以模擬基于全河的時間位移雙變量的月徑流二維矩陣，為流域水資源統(tǒng)一調(diào)度提供更為有力的依據(jù)。

為解決上述技術問題，本發(fā)明提供一種基于全河的月徑流時間位移二維矩陣的模擬方法，其特征是，包括如下步驟：

(1)構(gòu)建一種基于阿基米德copula及人工神經(jīng)網(wǎng)絡的三層模擬模型；

(2)構(gòu)建基于對稱阿基米德copula的單站點月徑流時間序列模擬模型以及基于非對稱阿基米德copula的多站點月徑流聯(lián)合分布模擬模型；

(3)用kendall系數(shù)采轉(zhuǎn)換法計算參數(shù)；

(4)構(gòu)建基于人工神經(jīng)網(wǎng)絡的沿河月徑流位移序列的模擬模型。

進一步的，所述步驟(1)構(gòu)建一種基于阿基米德copula及人工神經(jīng)網(wǎng)絡的三層模擬模型：

采用阿基米德copula模擬樣本徑流值，采用人工神經(jīng)網(wǎng)絡進行降階，求得模型參數(shù)與位移的關系函數(shù)，從而得到全河的二維徑流模擬函數(shù)：

式中fi，j(·；x)為徑流邊際分布，c(·；x)為copula函數(shù)，x為時間及位移二維變量，表示月徑流的徑流大小，i代表月徑流的時間，j代表月徑流的站點。

進一步的，根據(jù)研究區(qū)水文站的數(shù)量和研究資料的長度建立各站點各月的徑流變量矩陣代表月份；j＝1,2,…,n，n代表具體站點數(shù)目；表示yi,j變量的徑流值，k＝1,2,…,m，m代表時間長度。

進一步的，所述步驟(2)構(gòu)建基于對稱阿基米德copula的單站點月徑流時間序列模擬模型以及基于非對稱阿基米德copula的多站點月徑流聯(lián)合分布模擬模型：

月徑流自相關結(jié)構(gòu)沒有明顯的嵌套特點，假設當月徑流與前兩月徑流有相依關系，采用對稱阿基米德copula結(jié)構(gòu)，計算相鄰三個月的徑流聯(lián)合分布：

式中ui,j＝fi,j(·；x)；p表示概率。

根據(jù)下游站點徑流與其上游其余站點徑流均有相依關系這一相關特征，選取完全嵌套的非對稱阿基米德copula分層描述各站點月徑流與其上游所有站點月徑流聯(lián)合分布的相依結(jié)構(gòu)，其聯(lián)合分布為：

式中為不同層copula函數(shù)的發(fā)生器，n^*表示站點的序號，n^*小于等于n，為的逆函數(shù)。

進一步的，所述步驟(3)采用kendall系數(shù)轉(zhuǎn)換法計算參數(shù)：

通過變形，式(1.2)連續(xù)三個月的copula徑流聯(lián)合分布可表示為：

從而將三維copula函數(shù)轉(zhuǎn)換成二元copula函數(shù)計算，而完全嵌套的非對稱阿基米德copula是由一串二元copula嵌套而成，

因此，模型參數(shù)均可采用kendall系數(shù)轉(zhuǎn)換法計算參數(shù)，

copula函數(shù)參數(shù)與kendall系數(shù)的關系如下：

kendall系數(shù)計算方程如下：

式中τ¹ij用于時間層，τ²ij用于位移層，且：

進一步的，所述步驟(4)構(gòu)建基于人工神經(jīng)網(wǎng)絡的沿河月徑流位移序列的模擬模型：

采用aic法則分別對單站點月徑流時間序列模擬模型與多站點月徑流聯(lián)合分布模擬模型中不同阿基米德copula函數(shù)包括claytoncopula、gumbelcopula、frankcopula的模擬結(jié)果進行誤差分析，選出最佳copula函數(shù)類型：

式中v為參數(shù)數(shù)目，為均方差，w為樣本數(shù)量。

對上述選定的最佳copula函數(shù)的多站點模擬序列進行降階，將多站點月徑流聯(lián)合分布情況與沿河各點的關系函數(shù)轉(zhuǎn)化為模型中參數(shù)與各點距河源位移的關系函數(shù)：

ψχ＝gk(f²)(1.9)

式中f²為人工神經(jīng)網(wǎng)的輸出，gk(·)為轉(zhuǎn)換函數(shù)，ψχ為基于相應神經(jīng)元數(shù)目所求得的模型參數(shù)，χ為神經(jīng)元數(shù)目，

選擇一個隱含層的bp神經(jīng)網(wǎng)絡描述其關系，采用交叉檢驗的方法得到隱含層最佳神經(jīng)元的數(shù)目，確定神經(jīng)網(wǎng)絡結(jié)構(gòu)后求得沿河任意點相應對的參數(shù)，從而得到其徑流模擬值，

從而式(1.1)可表示為：

本發(fā)明所達到的有益效果：

豐富了中長期水文預報的模型庫，該模型不僅可模擬單站點的月徑流時間序列、多站點的月經(jīng)列聯(lián)合分布，還可以模擬時間位移雙變量的月徑流二維矩陣；模型中采用的對稱阿基米德copula函數(shù)和非阿基米德copula函數(shù)可以分別合理地描述時間和空間相關結(jié)構(gòu)；將所求沿河各點的聯(lián)合分布結(jié)構(gòu)轉(zhuǎn)化為計算分布參數(shù)與位移的關系，并通過bp神經(jīng)網(wǎng)絡構(gòu)造其關系函數(shù)。

附圖說明

圖1是本發(fā)明的流程圖；

圖2是對稱及非對稱copula函數(shù)結(jié)構(gòu)圖；

圖3是案例流域圖；

圖4是唐乃亥站點月徑流模擬結(jié)果圖；

圖5-1～5-12是α參數(shù)降尺度1～12月模擬結(jié)果圖；

圖6-1～6-12是β參數(shù)降尺度1～12月模擬結(jié)果圖；

圖7-1～7-12是θ參數(shù)降尺度1～12月模擬結(jié)果圖；

圖8是黃河干流月徑流時空二維模擬結(jié)果圖。

具體實施方式

下面結(jié)合附圖對本發(fā)明作進一步描述。以下實施例僅用于更加清楚地說明本發(fā)明的技術方案，而不能以此來限制本發(fā)明的保護范圍。

本發(fā)明是一種基于全河的月徑流時間位移二維矩陣的模擬方法，如圖所示1所示具體實施方式如下：

(1)構(gòu)建月徑流邊際分布擬合模型。

根據(jù)研究區(qū)水文站的數(shù)量和研究資料的長度建立各站點各月的徑流變量矩陣代表月份；j＝1,2,…,n，n代表具體站點數(shù)目；表示yi,j變量的徑流值，k＝1,2,…,m，m代表時間長度，選取pⅲ分布對各站點月徑流序列進行擬合，并采用極大似然法進行參數(shù)計算：

式中θ為參數(shù)向量。進而求得矩陣中各變量值。

(2)構(gòu)建單站點月徑流時間序列模擬模型。

氣候變化條件下河川徑流具有非線性、偏態(tài)分布等特點，可采用copula函數(shù)逐月構(gòu)造聯(lián)合分布函數(shù)，分別研究變量的邊緣分布和相關性結(jié)構(gòu)，不產(chǎn)生信息失真。根據(jù)sklar定理，令h(·,·)為具有邊緣分布f(·)和g(·)的聯(lián)合分布函數(shù)，且f(·)和g(·)連續(xù)，則存在唯一一個copula函數(shù)c(·,·)，滿足：

月徑流自相關結(jié)構(gòu)沒有明顯的嵌套特點，可以采用如圖2所示的對稱阿基米德copula結(jié)構(gòu)，計算相鄰三個月的徑流聯(lián)合分布，例如一月與二月與三月、十二月與次年一月與次年二月等。高維阿基米德copula只含有一個參數(shù)，但用計算過程復雜，可將其轉(zhuǎn)化為二元copula可簡化其參數(shù)求解：

則連續(xù)三個月的徑流聯(lián)合分布可表示為：

copula函數(shù)c(·,·)的具體函數(shù)類型多樣，其中阿基米德copula函數(shù)因具有完整的發(fā)生器而最為常用，其類型包括claytoncopula、gumbelcopula、frankcopula等。分別采用這三種二元阿基米德copula函數(shù)計算變量間的聯(lián)合分布從而推求連續(xù)三個月的徑流聯(lián)合分布。二元claytoncopula的聯(lián)合分布表達式為：

二元gumbelcopula的聯(lián)合分布表達式為：

二元frankcopula的聯(lián)合分布表達式為：

已知u1,u2,…,ui-1，可求得聯(lián)合分布的條件分布：

徑流矩陣中，隨機生成服從(0，1)均勻分布隨機數(shù)εt，并已知u1,1，可根據(jù)式(2.8)推求u2,1,u3,1，重復以上步驟順次推求其余月徑流。

(3)構(gòu)建基于多站點月徑流聯(lián)合分布模擬模型。

仍采用步驟(1)中各站點各月徑流邊際分布擬合結(jié)果，根據(jù)下游站點徑流與其上游其余站點徑流均有相依關系這一相關特征，選取如圖2所示的完全嵌套的非對稱阿基米德copula分層描述各站點月徑流與其上游所有站點月徑流聯(lián)合分布的相依結(jié)構(gòu)，例如站點1與站點2、站點3與站點1和站點2的聯(lián)合分布：

式中為不同層copula函數(shù)的發(fā)生器。

徑流矩陣中，隨機生成服從(0，1)均勻分布隨機數(shù)εt，并已知u1,1，可根據(jù)式(2.8)推求u1,2，根據(jù)式(2.5)、式(2.6)、式(2.7)可得u1,1,u1,2的聯(lián)合分布，繼而推求u1,3，順次推求其余站點月徑流。

(4)計算單站點月徑流時間序列模擬模型與多站點月徑流聯(lián)合分布模擬模型參數(shù)。

計算copula函數(shù)參數(shù)的方法有很多種，其中利用kendall系數(shù)求解參數(shù)的方法是最為簡單的，但該方法只適用于二元copula的情形。本模型本身結(jié)構(gòu)或結(jié)構(gòu)通過轉(zhuǎn)化即可符合kendall轉(zhuǎn)換法的使用條件。claytoncopula函數(shù)中參數(shù)與kendall系數(shù)的關系函數(shù)為：

gumbelcopula函數(shù)中參數(shù)與kendall系數(shù)的關系函數(shù)為：

frankcopula函數(shù)中參數(shù)與kendall系數(shù)的關系函數(shù)為：

數(shù)據(jù)樣本的kendall系數(shù)可按式(2.12)求得：

式中τ¹i,j用于時間層，τ²i,j用于位移層，且

(5)選取模擬單站點月徑流時間序列與多站點月徑流聯(lián)合分布的最佳copula類型

采用aic法則分別對單站點月徑流時間序列模擬模型與多站點月徑流聯(lián)合分布模擬模型中不同阿基米德copula函數(shù)的模擬結(jié)果進行誤差分析，選出最佳copula函數(shù)類型：

式中v為參數(shù)數(shù)目，為均方差，w為樣本數(shù)量。

(6)構(gòu)建沿河月徑流位移序列的模擬模型。

對上述選定的最佳copula函數(shù)的多站點模擬序列進行降階，將多站點月徑流聯(lián)合分布情況與沿河各點的關系函數(shù)轉(zhuǎn)化為模型中參數(shù)與各點距河源位移的關系函數(shù)。選擇一個隱含層的bp神經(jīng)網(wǎng)絡描述其關系，隱含層采用雙曲正切函數(shù)：

f¹l＝tansig(w¹lx+b¹l)(2.16)

輸出層采用線性傳遞函數(shù)：

采用rmse作為誤差函數(shù)：

(7)確定隱含層神經(jīng)元的個數(shù)。

采用交叉檢驗的方法得到隱含層最佳神經(jīng)元的數(shù)目，通過將樣本分成若干組，分別將每一組作為檢驗組，其他組作為訓練組，對不同神經(jīng)元數(shù)目下的結(jié)果進行檢驗，選擇誤差最小的情況時神經(jīng)元的數(shù)目。

之后確定神經(jīng)網(wǎng)絡結(jié)構(gòu)后求得沿河任意點相對應的參數(shù)，從而得到其徑流模擬值。

實施案例，

本實施案例選取黃河干流月徑流進行二維模擬。利用如圖3所示11個站點1956年至2010年的月徑流資料，該11個站點的基本資料如下表所示。

(1)構(gòu)建月徑流邊際分布擬合模型。

建立各站點各月徑流矩陣選取pⅲ分布對各站點月徑流序列進行擬合，并按式(2.1)進行參數(shù)計算。

(2)構(gòu)建單站點月徑流時間序列模擬模型。

依次按式(2.3)至式(2.8)計算一月與二月與三月、二月與三月與四月、三月與四月與五月、四月與五月與六月、五月與六月與七月、六月與七月與八月、七月與八月與九月、八月與九月與十月、九月與十月與十一月、十月與十一月與十二月、十二月與次年一月與次年二月連續(xù)三個月的徑流聯(lián)合分布。

隨機生成服從(0，1)均勻分布隨機數(shù)εt，并已知u11，可根據(jù)式(2.8)推求u21,u31，重復以上步驟順次推求其余月徑流。

(3)構(gòu)建基于多站點月徑流聯(lián)合分布模擬模型。

依次按式(2.9)計算各站點月徑流與其上游所有站點月徑流聯(lián)合分布的相依結(jié)構(gòu)，例如站點1與站點2、站點3與站點1和站點2的聯(lián)合分布。

隨機生成服從(0，1)均勻分布隨機數(shù)εt，并已知ui1，可根據(jù)式(2.8)推求ui2，根據(jù)式(2.5)、式(2.6)、式(2.7)可得ui1,ui2的聯(lián)合分布，繼而推求ui3，順次推求其余站點月徑流。

(4)計算單站點月徑流時間序列模擬模型與多站點月徑流聯(lián)合分布模擬模型參數(shù)。

按式(2.13)、式(2.14)分別計算相應數(shù)據(jù)間kendall系數(shù)，按式(2.10)至式(2.12)計算兩個模型的參數(shù)。將參數(shù)代入步驟三從而得到唐乃亥站點月徑流的模擬值如圖4所示。

(5)選取模擬單站點月徑流時間序列與多站點月徑流聯(lián)合分布的最佳copula類型

采用aic法則對單站點月徑流時間序列模擬模型中不同阿基米德copula函數(shù)的模擬結(jié)果進行誤差分析如下表所示，選出最佳copula函數(shù)類型為claytoncopula：

與多站點月徑流聯(lián)合分布模擬模型中不同阿基米德copula函數(shù)的模擬結(jié)果進行誤差分析如下表所示，選出最佳copula函數(shù)類型為gumbelcopula：

(6)構(gòu)建沿河月徑流位移序列的模擬模型。

對上述選定的最佳copula函數(shù)的多站點模擬序列進行降階，將多站點月徑流聯(lián)合分布情況與沿河各點的關系函數(shù)轉(zhuǎn)化為模型中參數(shù)與各點距河源位移的關系函數(shù)。選擇一個隱含層的bp神經(jīng)網(wǎng)絡描述其關系，隱含層采用雙曲正切函數(shù)，輸出層采用線性傳遞函數(shù)，采用交叉檢驗的方法得到隱含層最佳神經(jīng)元的數(shù)目為25個，之后求得沿河任意點相對應的參數(shù)如圖5、圖6、圖7所示，從而得到黃河干流的二維徑流模擬值如圖8所示。

以上所述僅是本發(fā)明的優(yōu)選實施方式，應當指出，對于本技術領域的普通技術人員來說，在不脫離本發(fā)明技術原理的前提下，還可以做出若干改進和變形，這些改進和變形也應視為本發(fā)明的保護范圍。