本發(fā)明屬于工程水文地質學中基流分割或水源劃分的技術領域。

背景技術：

基流分割是水文地質學中的一個重要問題，它在流域產(chǎn)匯流模型的模擬、流域水文分析計算、區(qū)域水資源評價和調(diào)查等領域具有十分重要的意義。這里的基流分割即將河川徑流的地表徑流、基流分割開。徑流是指在一定匯水區(qū)域上，由降水所形成的，在重力的作用下產(chǎn)生側向流動并沿著一定路徑匯集到河槽中的水流，其中，由地面匯集的水流稱為地表徑流，在飽和土層及巖石孔隙匯集的水流稱為地下徑流。水源劃分是基流分割的基礎，也是深入研究產(chǎn)匯流規(guī)律的核心性基礎性科學問題。水源組成及匯流過程示意圖見圖1，圖中，A為起漲點，B為地面徑流終止點，C為洪峰點，ACBD為河川徑流過程線，ACBA所包圍的面積為地面徑流量，ABDD′A′A所包圍的面積為基流。其中，N即為峰后地面徑流匯流時間。流域地表徑流經(jīng)驗公式最早是由Linsley等在《工程水文學》中提出的，其經(jīng)驗公式是：

N＝0.8F^0.2 (1)

式中，N——峰后地表徑流平均匯流時間，d；

F——流域面積(以平方公里記)，km²。

關于如何進行基流分割，前人提出了諸多方法，大致可分為圖解法、水文模型法、物理化學法、數(shù)學物理法及數(shù)值模擬法五類。圖解法是傳統(tǒng)的基流分割方法之一，其主觀性比較強且要求操作者熟悉流域的水文地質特征；水文模型法利用水文循環(huán)模型，基于質量守恒原理計算各水均衡項，物理意義明確，具有客觀可重復性的特點，但對使用者能力具有較高的要求，且模型所需調(diào)試的參數(shù)較多，模型參數(shù)的選取困難。物理化學法利用同位素放射性示蹤劑進行示蹤，結果準確，可信度高，但對監(jiān)測儀器和監(jiān)測區(qū)水文地質條件要求比較高且所需費用較高。數(shù)學物理法通過建立數(shù)學物理方程，較好刻畫了基流產(chǎn)流的物理機制，但當結果出現(xiàn)較大偏差時，多參數(shù)調(diào)整與適用范圍確定的難度較大。數(shù)值模擬法基于信號處理技術，通過河道徑流量計算基流量，其中，主要方法包括滑動最小值法、濾波法和加里寧系列方法，具有精度高，客觀性強，操作簡便，計算速度快，可重復使用等特點，并且可直接易于采用計算機自動實現(xiàn)，在實踐中得到了廣泛應用和認可。目前，在進行基流分割時，數(shù)值模擬法應用最廣，其中數(shù)字濾波法是近年來國際上研究最多的基流分割方法，在美國65個流域應用Eckhardt數(shù)字濾波法與其他方法進行比較，研究表明該方法進行地下徑流分割最為合理?？紤]到Eckhardt遞歸數(shù)字濾波方程中參數(shù)I_GF.max取值僅僅根據(jù)流域的下墊面條件和水文條件簡單地取0.80、0.50和0.25難免會產(chǎn)生較大的誤差，因此在取值的問題上做了一些改進，采用滑動最小值法進行估算對Eckhardt遞歸數(shù)字濾波方程中參數(shù)I_GF.max值進行確定，這樣可以避免參數(shù)取值的隨意性，提高計算結果的精確度。但是在滑動最小值估算過程中要試算流域峰后地面徑流匯流時間N，其試算方法人為主觀性比較強，誤差較大，容易對結果造成影響。

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

本發(fā)明提供一種針對半濕潤半干旱地區(qū)及濕潤地區(qū)的基流分割方法，以解決水文地質學中基流分割工程問題。

本發(fā)明采取的技術方案是：包括下列步驟：

步驟1，收集數(shù)據(jù)，包括日徑流量Q、日降雨量P和匯水區(qū)流域面積F，以平方公里記，km²；

根據(jù)日徑流量資料，采用退水曲線法將標準退水曲線移動到透明紙上，再將其覆蓋到要分割的流量過程線的退水段上，注意比例尺一致，基線與時間軸重合，左右平移使兩者退水段尾部吻合，則兩線開始重疊的時刻，就可以作為地面徑流的終止點，地面徑流終止點距離實測流量過程線的洪峰流量出現(xiàn)的時刻的時距為N,即得到峰后地表徑流平均匯流時間，天數(shù)d；

以退水曲線求得的N值和匯水區(qū)流域面積F作為初始樣本，流域地表徑流經(jīng)驗公式如下：

N＝0.8F^0.2 (1)

步驟2，流域峰后地面徑流匯流時間經(jīng)驗公式率定

采用遺傳算法對經(jīng)驗公式(1)中的參數(shù)進行優(yōu)化，完成優(yōu)化經(jīng)驗公式的率定，此過程通過MATLAB運行遺傳算法參數(shù)優(yōu)化程序完成，具體步驟如下：

(a)染色體編碼

為了簡明起見，將公式(1)作如下轉化：

這里的x₁和x₂為經(jīng)驗公式的參數(shù)，在遺傳算法中作為染色體，F(xiàn)為經(jīng)驗公式的自變量，即流域面積，在遺傳算法中作為輸入，變量x₁和x₂的取值空間取[0，1]，染色體采用浮點數(shù)編碼；

(b)目標函數(shù)

根據(jù)經(jīng)驗公式推求的基本原則，遺傳算法的目標函數(shù)可以表示為：

式中，N'——經(jīng)驗公式N(x₁,x₂)推求的流域地面徑流平均匯流時間，d；

N——流域地面徑流平均匯流時間計算值，d；

對公式(3)搜尋某一組參數(shù)使經(jīng)驗公式推求的流域地面徑流平均匯流時間N'與計算值N差的平方倒數(shù)最大，為了防止目標函數(shù)的計算值過小，引入擴大系數(shù)M；公式(3)轉換為：

(c)運行過程

遺傳算法優(yōu)化經(jīng)驗公式參數(shù)的過程從一組隨機產(chǎn)生的初始解，即“種群”，開始進行搜索，種群中的每一個個體，即問題的一個解，稱為“染色體”；遺傳算法通過染色體的“適應值”來評價染色體的好壞，適應值大的染色體被選擇的幾率高，相反，適應值小的染色體被選擇的可能性小，被選擇的染色體進入下一代；下一代中的染色體通過交叉和變異等遺傳操作，產(chǎn)生新的染色體，即“后代”；經(jīng)過若干后代之后，算法收斂于最好的染色體，該染色體就是問題的最優(yōu)解或近似解；

(d)遺傳算法優(yōu)化模型參數(shù)

遺傳算法參數(shù)有一個取值范圍，在其范圍內(nèi)確定其參數(shù)值如表1；

表1 遺傳算法參數(shù)

經(jīng)過遺傳算法優(yōu)化參數(shù)后得到優(yōu)化經(jīng)驗公式為：

N＝0.1200F^0.3055 (5)

步驟3，流域地面匯流時間優(yōu)化經(jīng)驗公式的改進。

對于半干旱半濕潤地區(qū)、及濕潤地區(qū)，雨強仍然是制約匯流時間的一個重要指標，進而，對原經(jīng)驗公式(1)加以改進，引入變量平均雨強P，建立流域地面徑流平均匯流時間與流域面積與平均雨強之間的相關關系，

再采用遺傳算法進行參數(shù)優(yōu)化，具體步驟如下：

(a)染色體編碼

這里的x₁、x₂和x₃為經(jīng)驗公式的參數(shù)，在遺傳算法中作為染色體，F(xiàn)、P為經(jīng)驗公式的自變量，即流域面積和日平均雨強，在遺傳算法中作為輸入，變量x₁、x₂和x₃的取值空間取[0，1]，染色體采用浮點數(shù)編碼；

(b)目標函數(shù)

根據(jù)經(jīng)驗公式推求的基本原則，遺傳算法的目標函數(shù)可以表示為：

式中，N″——經(jīng)驗公式N(x₁,x₂,x₃)推求的流域地面徑流平均匯流時間，d；

N——流域地面徑流平均匯流時間計算值，d；

對公式(8)搜尋某一組參數(shù)使經(jīng)驗公式推求的流域地面徑流平均匯流時間N″與計算值N差的平方倒數(shù)最大，為了防止目標函數(shù)的計算值過小，引入擴大系數(shù)M′，公式(8)轉換為：

(c)運行過程

遺傳算法優(yōu)化經(jīng)驗公式參數(shù)的過程從一組隨機產(chǎn)生的初始解，即“種群”，開始進行搜索，種群中的每一個個體，即問題的一個解，稱為“染色體”；遺傳算法通過染色體的“適應值”來評價染色體的好壞，適應值大的染色體被選擇的幾率高，相反，適應值小的染色體被選擇的可能性小，被選擇的染色體進入下一代；下一代中的染色體通過交叉和變異等遺傳操作，產(chǎn)生新的染色體，即“后代”；經(jīng)過若干后代之后，算法收斂于最好的染色體，該染色體就是問題的最優(yōu)解或近似解；

(d)遺傳算法優(yōu)化模型參數(shù)

一般情況下，遺傳算法參數(shù)有一個取值范圍，我們在其范圍內(nèi)確定其參數(shù)值如表2；

表2 遺傳算法參數(shù)

經(jīng)過遺傳算法優(yōu)化參數(shù)后得到改進的優(yōu)化經(jīng)驗公式為：

N＝0.0080F^0.3206P^1.4824 (10)

步驟4，I_GF.max值的確定

首先采用滑動最小值法計算出N為1-25所對應的每一年的地下徑流指數(shù)I_GF，此過程通過BFI源程序完成，根據(jù)改進的經(jīng)驗公式(10)計算N值，確定其N值后，找到N值所對應的每一年的地下徑流指數(shù)I_GF，選取這些年中最大的地下徑流指數(shù)I_GF，即為I_GF.max；

步驟5，基流分割

采用Eckhardt遞歸數(shù)字濾波法對徑流進行分割，Eckhardt遞歸數(shù)字濾波方程：

式中a為濾波參數(shù),I_GF.max為最大地下徑流指數(shù)(groundwaterflow index)，表示河川徑流中地下徑流所占比例，

這里參數(shù)a取0.98，I_GF.max即為上述步驟4所得。

q_t與q_t-1分別表示第t和t-1時刻的地面徑流，以日為時間步長，q_t≥0，m³/s；

Q_t分別表示第t時刻的總徑流，即河川徑流，m³/s。

本發(fā)明對于半濕潤半干旱地區(qū)及濕潤地區(qū)和濕潤地區(qū)而言，引入變量雨強，改進了經(jīng)驗公式，提高了計算精度，針對半濕潤半干旱地區(qū)及濕潤地區(qū)和濕潤地區(qū)提出了一種基流分割方法。具體采用濾波法來進行基流分割，其中，峰后地面徑流匯流時間參數(shù)采用改進經(jīng)驗公式推求，并采用遺傳算法來優(yōu)化經(jīng)驗公式參數(shù)，參數(shù)選擇具有全局最優(yōu)性。

本發(fā)明的優(yōu)點是在使用Eckhardt遞歸數(shù)字濾波法進行基流分割中，需要用滑動最小值法對Eckhardt遞歸數(shù)字濾波方程中參數(shù)I_GF.max值進行確定，而滑動最小值法估算過程需要計算峰后地面徑流匯流時間N，之前是采用試算法進行估算，但其試算方法人為主觀性比較強，誤差較大，容易對結果造成影響。本發(fā)明為了克服滑動最小值試算的缺點，針對選定的流域得到一個流域峰后地面徑流匯流時間經(jīng)驗公式，這樣就能直接計算峰后地面徑流匯流時間N，這種方法針對選定的流域有一個固定的經(jīng)驗公式，不僅簡便，而且又可以提高精確度，其中，經(jīng)驗公式又有以下兩點改進：

1)流域峰后地面徑流匯流時間經(jīng)驗公式率定，Linsley定義原始經(jīng)驗公式時，只是粗略觀察得到地表徑流匯流時間N，比較粗糙，且精度較低，本發(fā)明為了改進這個缺陷，采用遺傳算法去優(yōu)化原始經(jīng)驗公式，使其能夠更好的使用于實際工作中。

2)流域地面匯流時間經(jīng)驗公式改進，原始的經(jīng)驗公式只是考慮到面積對匯流時間的影響，并沒有考慮到雨強，而對于半濕潤半干旱地區(qū)及濕潤地區(qū)而言，雨強仍然是制約匯流時間的一個重要指標，本發(fā)明引入變量雨強，建立了新的經(jīng)驗公式，提高了計算精度。

附圖說明

圖1是水源組成示意圖；

圖2是步驟2)中遺傳算法優(yōu)化經(jīng)驗公式參數(shù)的過程圖；

圖3是步驟3)中遺傳算法優(yōu)化經(jīng)驗公式參數(shù)的過程圖；

圖4是嫩江流域計算結果比較圖；

圖5是嫩江流域計算結果比較圖。

具體實施方式

包括下列步驟：

步驟1，收集數(shù)據(jù)，包括日徑流量Q、日降雨量P和匯水區(qū)流域面積F，以平方公里記，km²；

根據(jù)日徑流量資料，采用退水曲線法將標準退水曲線移動到透明紙上，再將其覆蓋到要分割的流量過程線的退水段上，注意比例尺一致，基線與時間軸重合，左右平移使兩者退水段尾部吻合，則兩線開始重疊的時刻，就可以作為地面徑流的終止點，地面徑流終止點距離實測流量過程線的洪峰流量出現(xiàn)的時刻的時距為N,即得到峰后地表徑流平均匯流時間，天數(shù)d；參見圖1；

以退水曲線求得的N值和匯水區(qū)流域面積F作為初始樣本，流域地表徑流經(jīng)驗公式如下：

N＝0.8F^0.2 (1)

步驟2，流域峰后地面徑流匯流時間經(jīng)驗公式率定

采用遺傳算法對經(jīng)驗公式(1)中的參數(shù)進行優(yōu)化，完成優(yōu)化經(jīng)驗公式的率定，此過程通過MATLAB運行遺傳算法參數(shù)優(yōu)化程序完成，具體步驟如下：

(a)染色體編碼

為了簡明起見，將公式(1)作如下轉化：

這里的x₁和x₂為經(jīng)驗公式的參數(shù)，在遺傳算法中作為染色體，F(xiàn)為經(jīng)驗公式的自變量，即流域面積，在遺傳算法中作為輸入，變量x₁和x₂的取值空間取[0，1]，染色體采用浮點數(shù)編碼；

(b)目標函數(shù)

根據(jù)經(jīng)驗公式推求的基本原則，遺傳算法的目標函數(shù)可以表示為：

式中，N′——經(jīng)驗公式N(x₁,x₂)推求的流域地面徑流平均匯流時間，d；

N——流域地面徑流平均匯流時間計算值，d；

對公式(3)搜尋某一組參數(shù)使經(jīng)驗公式推求的流域地面徑流平均匯流時間N′與計算值N差的平方倒數(shù)最大，為了防止目標函數(shù)的計算值過小，引入擴大系數(shù)M；公式(3)轉換為：

(c)運行過程

遺傳算法優(yōu)化經(jīng)驗公式參數(shù)的過程可以采用圖2描述，從一組隨機產(chǎn)生的初始解，即“種群”，開始進行搜索，種群中的每一個個體，即問題的一個解，稱為“染色體”；遺傳算法通過染色體的“適應值”來評價染色體的好壞，適應值大的染色體被選擇的幾率高，相反，適應值小的染色體被選擇的可能性小，被選擇的染色體進入下一代；下一代中的染色體通過交叉和變異等遺傳操作，產(chǎn)生新的染色體，即“后代”；經(jīng)過若干后代之后，算法收斂于最好的染色體，該染色體就是問題的最優(yōu)解或近似解；

(d)遺傳算法優(yōu)化模型參數(shù)

遺傳算法參數(shù)有一個取值范圍，在其范圍內(nèi)確定其參數(shù)值如表1；

表1 遺傳算法參數(shù)

經(jīng)過遺傳算法優(yōu)化參數(shù)后得到優(yōu)化經(jīng)驗公式為：

N＝0.1200F^0.3055 (5)

步驟3，流域地面匯流時間經(jīng)驗公式改進。

對于半干旱半濕潤地區(qū)、及濕潤地區(qū)，雨強仍然是制約匯流時間的一個重要指標，進而，對原經(jīng)驗公式(1)加以改進，引入變量平均雨強P，建立流域地面徑流平均匯流時間與流域面積與平均雨強之間的相關關系，

再采用遺傳算法進行參數(shù)優(yōu)化，具體步驟如下：

(a)染色體編碼

這里的x₁、x₂和x₃為經(jīng)驗公式的參數(shù)，在遺傳算法中作為染色體，F(xiàn)、P為經(jīng)驗公式的自變量，即流域面積和日平均雨強，在遺傳算法中作為輸入，變量x₁、x₂和x₃的取值空間取[0，1]，染色體采用浮點數(shù)編碼；

(b)目標函數(shù)

根據(jù)經(jīng)驗公式推求的基本原則，遺傳算法的目標函數(shù)可以表示為：

式中，N″——經(jīng)驗公式N(x₁,x₂,x₃)推求的流域地面徑流平均匯流時間，d；

N——流域地面徑流平均匯流時間計算值，d；

對公式(8)搜尋某一組參數(shù)使經(jīng)驗公式推求的流域地面徑流平均匯流時間N″與計算值N差的平方倒數(shù)最大，為了防止目標函數(shù)的計算值過小，引入擴大系數(shù)M′，公式(8)轉換為：

(c)運行過程

遺傳算法優(yōu)化經(jīng)驗公式參數(shù)的過程可以采用圖3描述，從一組隨機產(chǎn)生的初始解，即“種群”，開始進行搜索，種群中的每一個個體，即問題的一個解，稱為“染色體”；遺傳算法通過染色體的“適應值”來評價染色體的好壞，適應值大的染色體被選擇的幾率高，相反，適應值小的染色體被選擇的可能性小，被選擇的染色體進入下一代；下一代中的染色體通過交叉和變異等遺傳操作，產(chǎn)生新的染色體，即“后代”；經(jīng)過若干后代之后，算法收斂于最好的染色體，該染色體就是問題的最優(yōu)解或近似解；

(d)遺傳算法優(yōu)化模型參數(shù)

一般情況下，遺傳算法參數(shù)有一個取值范圍，我們在其范圍內(nèi)確定其參數(shù)值如表2；

表2 遺傳算法參數(shù)

經(jīng)過遺傳算法優(yōu)化參數(shù)后得到優(yōu)化經(jīng)驗公式為：

N＝0.0080F^0.3206P^1.4824 (10)

步驟4，I_GF.max值的確定

首先采用滑動最小值法計算出N為1-25所對應的每一年的地下徑流指數(shù)I_GF，此過程通過BFI源程序完成，根據(jù)改進的經(jīng)驗公式(10)計算N值，確定其N值后，找到N值所對應的每一年的地下徑流指數(shù)I_GF，選取這些年中最大的地下徑流指數(shù)I_GF，即為I_GF.max；

步驟5，基流分割

采用Eckhardt遞歸數(shù)字濾波法對徑流進行分割，Eckhardt遞歸數(shù)字濾波方程：

式中a為濾波參數(shù),I_GF.max為最大地下徑流指數(shù)(groundwaterflow index)，表示河川徑流中地下徑流所占比例，

這里參數(shù)a取0.98，I_GF.max即為上述步驟4所得；

q_t與q_t-1分別表示第t和t-1時刻的地面徑流，以日為時間步長，q_t≥0，m³/s；

Q_t分別表示第t時刻的總徑流，即河川徑流，m³/s；

此過程通過計算網(wǎng)址完成，其網(wǎng)址為：

https://engineering.purdue.edu/～what/。

檢驗例對流域地面匯流時間經(jīng)驗公式、優(yōu)化經(jīng)驗公式和改進的優(yōu)化經(jīng)驗公式的檢驗。

由于嫩江流域是半干旱半濕潤地區(qū)，所以選用嫩江流域5個匯水區(qū)作為作為檢驗區(qū)域。采用優(yōu)化經(jīng)驗公式(5)和改進的優(yōu)化經(jīng)驗公式(10)進行計算，估算結果見圖4和圖5，分別進行誤差分析，即分別用經(jīng)驗公式(1)、優(yōu)化經(jīng)驗公式(5)、改進的優(yōu)化經(jīng)驗公式(10)計算出來的結果和初始值N比較,其差值結果見表3。

表3 誤差統(tǒng)計

從表3中可以看到，引入變量雨強并優(yōu)化之后的改進優(yōu)化經(jīng)驗公式精度最高，優(yōu)化經(jīng)驗公式精度次之，原經(jīng)驗公式精度最低，最終采用公式(10)改進的優(yōu)化經(jīng)驗公式作為應用公式。

下邊通過本發(fā)明方法在松花江的應用來進一步說明本發(fā)明。

以白山匯水區(qū)為例，具體步驟如下：

1)收集數(shù)據(jù)，白山匯水區(qū)面積分別為19000km²和日平均雨強6.59mm；

2)采用公式(10)改進的優(yōu)化經(jīng)驗公式，求得流域地表徑流平均匯流時間，即：

N＝0.0080F^0.3206P^1.4824

＝0.0080×19000^0.3206×6.59^1.4824

＝3.0816

由此看出，N取3比較合理。

3)I_GF.max值的確定。

首先采用滑動最小值法進行地下徑流分割，分別計算出每一年的地下徑流指數(shù)I_GF，當N取3時，各年I_GF計算值見表4，由表4可知1971-1980年10a中1978年的I_GF最大，由此得I_GF.max為0.744

表4 I_GF值計算表

4)基流分割

Eckhardt遞歸數(shù)字濾波方程：

式中a為濾波參數(shù),I_GF.max為最大地下徑流指數(shù)(groundwaterflow index),表示河川徑流中地下徑流所占比例。

q_t為時段(1d)內(nèi)過濾后的地面徑流(快速響應)；

Q_t為原始流量；

采用Eckhardt遞歸數(shù)字濾波法對1971-1980年的徑流資料進行分析計算，這里a取0.98，I_GF.max取0.744，其地下徑流(基流)分割計算成果見表5和表6。

表5 第二松花江白山站1971-1980年月平均徑流量

表6 第二松花江白山站各年徑流分割成果表