月徑流隨機模擬序列優選方法

劉欣欣，馬細霞,2，程旭，王倩麗，張靜文

(1.鄭州大學水利科學與工程學院，鄭州 450001;2.鄭州大學黃河生態保護與區域協調發展研究院，鄭州 450001)

徑流序列作為水資源系統的輸入項，通常用來分析系統未來的發展趨勢。通過水文隨機模型生成足夠長的徑流序列，可以彌補實測資料代表性不足的缺陷，已廣泛應用于水庫(群)優化調度[1-2]、水資源系統風險分析[3-4]等領域。

目前，水文隨機模型的精度評價常采用“統計一致性”準則，即生成的模擬序列必須和實測資料具有相近的統計特征。如:Weiss等[5]以均值、均方差、滯后1和滯后2的自相關系數作為檢驗指標，對實測序列和模擬長序列的高分位數進行比較，評估季節性自回歸模型[SAR(p)]保持實測序列統計特征方面的性能。Hao等[6]以科羅拉多河為研究對象，通過分析月徑流模擬序列基本統計參數的相對誤差，發現平均值、標準差和偏度的相對誤差均小于10%，滯后1的自相關系數相對誤差均小于20%，由此得出了模擬序列能較好地保持實測序列截口統計特征的結論。Hamid等[7]采用平均絕對百分比誤差(EMAP)等方法來分析模擬結果的均值、標準差、偏態系數和最大值，根據EMAP的最小值表示最佳模式，優選出最適合的季節性自回歸綜合移動平均模型。高瑞忠等[8]采用平均相對誤差和平均絕對誤差方法，評價巴拉格爾河流域年季尺度下的徑流量模擬效果，據此得出了神經網絡方法對該流域徑流量進行模擬總體適用性較好的結論。傳統評價方法在保證模擬系列均值、方差等基本統計特征較為一致的情況下，能夠模擬出具有不同時空分配特性的徑流過程，對水庫(群)優化調度、水資源系統風險分析、水資源工程規模確定具有重要意義。

徑流過程是十分復雜的隨機過程，但也具有一定的統計規律，譬如年內分配過程具有汛期和非汛期的交替變化特征，同一地區汛期最大4個月所處月份基本固定，多年平均最大4個月水量占全年平均水量的比例也較為穩定。徑流模擬序列具有不同的時空分配特性，若僅保證模擬系列均值、方差等基本統計特征較為一致，則有可能存在某些年份的年內分配過程不合理問題，因此基于此模擬序列獲得的應用成果存在一定的安全風險。為此，馬細霞等[9]引入樣本熵指標對徑流序列的復雜特性進行檢驗，結果表明樣本熵指標能大體診斷出模擬序列內較大月徑流量和位置的異常情況，但其合理性沒有得到充分論證。

以往評價方法的不足之處在于未能全面檢驗模擬序列年內復雜變化的特性，因此，有必要增加一些反映徑流時序變化的特征指標，如連續最大4個月徑流百分率[10]、徑流年內分配不均勻系數、徑流年內集中度[10-11]等，在保證模擬序列年內分配多樣性的同時，排除年內分配異常的情況。此外，同一隨機模型所生成的模擬序列在時程變化上差異較大[9]，以往的評價方法在工程規劃設計和風險分析方面，未能篩選出最適宜的徑流模擬序列。鑒于此，本文提出截口-復雜性指標體系，以淮河流域河南段干支流上的典型水文研究站月徑流隨機模擬為應用實例，對季節性自回歸模型模擬序列的截口特性和復雜特性進行全面檢驗，并基于灰色關聯分析法對以往指標體系和截口-復雜性指標體系的計算結果進行優選和分析，以期為月徑流隨機模型模擬序列優選提供理論依據。

1 模擬序列評價指標體系及指標量化

1.1 評價指標體系構建

為全面評價模擬序列，本文從以下2個方面構建截口-復雜性指標體系。截口特征方面：均值、變差系數、偏差系數反映序列截口統計特性，相鄰截口滯時1、2的自相關系數反映徑流序列的時間相依特性。這些指標在隨機模擬序列檢驗中已得到廣泛利用，故本文仍將其作為月徑流隨機模擬效果評價的主要指標。復雜特性方面：河川徑流的時空分布規律是水資源評價的重要內容，也是河流類型和水利規劃的一項重要指標[12]，本文引入連續最大4個月徑流百分率、集中度和不均勻系數作為反映徑流變化過程的指標，以期更好地反映實測序列與模擬序列年內分配的復雜特征。樣本熵表示非線性動力學系統產生新信息的概率，主要用來定量地刻畫系統的規則度及復雜度[13]，當用于度量時間序列復雜性時，其值越大，序列產生新信息的概率越高，序列越復雜[14]。由于其較少地依賴于時間序列長度，所需數據序列短、計算成果更為穩定，已在水文氣象研究領域廣泛應用，因此可以作為描述徑流序列復雜性特征的指標[15]。傳統評價指標與截口-復雜性指標體系對比見表1。

1.2 評價指標量化方法

截口參數采用統計學方法計算[16]，以下主要介紹樣本熵、集中度和不均勻系數的計算方法。

1.2.1樣本熵計算

設原始數據為X(1),X(2),…,X(N)，共N點，可按如下步驟計算其樣本熵[17]。

按序號連續順序，組成一組m維矢量，從Xm(1)到Xm(N-m+1)，其中

Xm(i)=[x(i),x(i+1),…,x(i+m-1)],i=1～N-m+1

(1)

定義X(i)與X(j)之間的距離d[Xm(i),Xm(j)]為兩者對應元素中差值最大的一個，即

d[Xm(i),Xm(i)]=max(|x(i+k)-x(j+k)|)，

k=0～m-1;i,j-1～N-m+1;i≠j

(2)

給定閾值r，對每一個1≤i≤N-m值，統計d[Xm(i),Xm(j)]

Bm,i(r)={d[Xm(i),Xm(j)]

(3)

求其對所有i的平均值，即

(4)

將維數增加1，變為m+1維矢量，重復式(1)～(3)的步驟后，得

(5)

式中：Bm(r)和Bm+1(r)分別為m點和m+1點兩序列相似的概率，則序列理論上的樣本熵

(6)

當N為有限時，得出序列樣本熵估計值

H(m,r,N)=-ln[Bm+1(r)/Bm(r)]

(7)

參數m、r的選擇是樣本熵估計的關鍵，但目前尚無最佳標準，通常m=2，r=0.2SD(SD為原始序列的標準差)[17]。

1.2.2集中度計算

徑流在年內的集中程度用集中度(Cd)表示，計算公式[18]為

(8)

式中：R為多年平均月徑流計算得出的年徑流量；Rx、Ry分別為12個月份 (i=1,2,…,12)的分量之和所構成的水平、垂直分量，其中

(9)

式中：ri為第i月平均月徑流；θi為第i月份對應的方位角，以各月月中值代表的角度數值來表示。

1.2.3不均勻系數計算

不均勻系數(Ct)的計算公式[19]為

(10)

式中：Ct(t)為第t年的徑流年內分配不均勻系數；σ(t)為第t年的年內時段平均流量的均方差；Q(t)為第t年的平均流量，m3/s；λ為月份，λ=1,2,…,m，此處m=12。徑流年內分配不均勻系數反映了對徑流調控的難易程度。對于n年系列的不均勻系數的均方差，其值越大，則表示徑流不均勻性的年際變化越大。

1.3 評價指標分析方法

平均絕對百分比誤差(EMAP)屬于無量綱統計量，可以在不同變量間進行比較，其值越小，說明模型模擬或預測的精確度越高，因此可表征各評價指標中模擬序列的精度。EMAP的計算公式為

(11)

式中：Ai為實測序列第i個指標計算值；Fi為模擬序列第i個指標計算值。

2 隨機模擬序列的灰色系統評價方法

隨機模擬生成徑流序列最基本的假定是未來事件應具有和實測序列相同的隨機性質，這就要求生成的模擬序列必須和實測資料具有相近的統計特征，這也是優選隨機模型和評價生成序列的主要依據。灰色關聯分析法是一種衡量因素的發展態勢和特征變化相似程度的方法，通過灰色關聯度表征數據序列之間聯系的緊密程度。本文采用灰色關聯分析法，計算模擬序列與實測資料在截口、復雜性指標的綜合關聯度，優選較為理想的隨機模擬序列。灰色系統關聯分析的具體步驟如下[20]。

確定參考數列與比較數列。本文采用截口-復雜性指標體系中的9個指標作為月徑流序列的特征指標，由實測序列統計所得的9個特征指標值構成參考數列X0，第i個隨機模擬序列統計所得的9個特征指標值構成比較數列Xi(i=1,2,…,n)。

X0={x0(1),x0(2),…,x0(9)}

(12)

Xi={xi(1),xi(2),…,xi(9)}

(13)

式中：x0(1),x0(2),…,x0(9)分別為參考數列X0的9個特征指標值；xi(1),xi(2),…,xi(9)分別為比較數列Xi的9個特征指標值。

指標值的無量綱化處理。計算公式為

(14)

(15)

式中：x0(k)為參考數列X0的第k個特征指標值；xi(k)為比較數列Xi的第k個特征指標值；yi(k)為比較數列Xi第k個特征指標的效果測度；k=1,2,…,9。

求參考數列與比較數列的灰色關聯系數。計算公式為

ξi(k)=

(16)

式中：ξi(k)為參考數列X0與比較數列Xi第k個特征指標的關聯系數；ρ為分辨系數，表示平均分辨水平，其值為0～1，一般取0.5[20]。

求參考數列與比較數列的加權關聯度。為便于從整體上比較模擬序列，有必要將各指標的關聯系數集中為一個值，即加權關聯度為ri，其表達式為

(17)

式中：ωi(k)為比較數列Xi第k個特征指標的權重，考慮到截口特性、復雜變化特性方面的各指標重要性相當，故本文采用等權處理，即：ωi(k)=1/9。加權關聯度越大，說明相應的模擬序列在各指標的綜合評價中表現越好。

3 實例分析與討論

3.1 研究區概況及數據來源

淮河是河南省內流域面積最大的水系，境內控制流域面積8.83萬km2，占淮河流域總面積的32%，占河南省總面積的53%。為驗證本文方法的合理性和有效性，在淮河干流、沙河和史河支流上，分別選取流域面積相差較大的息縣(10 190 km2)、漯河(12 150 km2)和蔣家集水文研究站(5 930 km2)月徑流隨機模擬為應用實例，以期為淮河流域中上游徑流模擬序列選取提供理論支撐。以上3站實測年月徑流資料(1951—2018，共68年)均來源于《河南省水情手冊》，這些資料均通過了“三性”審查，可以用于模擬計算。

3.2 季節性自回歸模型建立

根據RIC定階方法[21]，計算漯河站、息縣站、蔣家集站SAR(p)模型不同月份的RIC值。其中：漯河站和息縣站12個月的RIC最小值對應階數均為1階；蔣家集站除11月份外，其余11個月RIC最小值對應的階數均為1階。因此，根據RIC 準則，3個水文站均采用SAR(1)模型進行隨機模擬。

3.3 隨機模擬序列精度分析

通過所建SAR(1)模型，隨機生成10組樣本容量為680年的模擬序列，樣本的選取具有隨機性。為便于分析，將漯河站10組模擬序列記為1-1、1-2……1-10，息縣站記為2-1、2-2……2-10，蔣家集站記為3-1、3-2……3-10。采用相對誤差和平均絕對百分比誤差(EMAP)，從截口、復雜變化方面分析模擬序列的擬合效果，并結合模擬序列多年平均及各水平年年內分配情況，驗證截口、復雜變化特性指標評價的合理性。

3.3.1截口統計特性精度分析

模擬序列截口指標的相對誤差和EMAP結果見表2。由表2可知，各水文站10組模擬序列在截口特性方面的精度存在差異，3個水文研究站的10組模擬序列的均值和變差系數相對誤差均在15%范圍內，偏態系數的誤差相對較大，一階自相關系數的誤差平均值在8%以內，二階自相關系數的誤差平均值在15%范圍內。EMAP計算結果表明:漯河站EMAP值最小的是1-5序列，其次是1-8、1-10、1-9序列；息縣站EMAP值最小的是2-1、2-10序列，其次是2-9序列；蔣家集站EMAP值最小的是3-4、3-7序列，其次是3-1序列。與同一水文研究站的其他模擬序列相比，這些模擬序列在截口統計特性評價中表現較好。

表2 模擬序列截口指標相對誤差及EMAP結果統計Tab.2 Statistics of RE and EMAP results of simulated sequence section index %

3.3.2復雜變化特性精度分析

模擬序列復雜變化指標的相對誤差和EMAP結果見表3。可以看出，同一水文研究站各模擬序列在復雜變化特性方面的精度存在差異。其中：3個水文站連續最大4個月徑流百分率、不均勻系數的相對誤差在13%范圍內；集中度、樣本熵指標變化幅度較大，說明模擬序列在集中程度和復雜特性方面差異較大。EMAP計算結果表明：漯河站EMAP值最小的是1-5序列，其次是1-2序列；息縣站EMAP值最小的是2-1序列，2-2、2-3序列次之；蔣家集站EMAP值最小的是3-7序列，其次是3-9序列。與同一水文研究站的其他模擬序列相比，這些模擬序列較好保持了實測序列在形態變化方面的特性。

表3 模擬序列復雜變化指標相對誤差及EMAP統計Tab.3 Statistics of RE and EMAP results of complex change index of simulation sequence %

3.3.3模擬序列精度綜合分析

為驗證截口-復雜性指標體系的合理性，采用距平百分率法對徑流序列進行豐、平、枯劃分，統計各站模擬序列多年平均及豐、平、枯水年年內分配的相對誤差，結果見圖1。可以發現：漯河站模擬序列豐、平水年以及多年平均月徑流的擬合效果相對較好，相對誤差均在20%范圍內，枯水年各序列間相對誤差差異較大，年內分配總體擬合較好的是1-5、1-7序列；息縣站模擬序列各水平年及多年平均的相對誤差在25%范圍內，豐水年的相對誤差均在12%范圍內，年內分配總體擬合較好的是2-1序列；蔣家集站模擬序列多年平均及各水平年的相對誤差整體在25%范圍內，各序列相對誤差差異明顯，年內分配總體擬合較好的是3-4、3-5、3-7序列。

圖1 各站模擬序列多年平均以及豐/平/枯水年相對誤差Fig.1 The relative errors of the simulated series of each station in annual average,wet,normal and dry years

對比模擬序列截口、復雜性精度分析結果與水平年精度分析結果 (見表4)，可以看出：同一水文研究站的模擬序列在截口、復雜變化特性方面擬合效果相差各異。例如：漯河站模擬序列1-8、1-9、1-10在截口統計特性評價中EMAP值較小，說明其較好的保持了實測序列截口統計特性，但在復雜變化特性和年內分配過程方面EMAP值均較大；蔣家集站3-7、3-9序列復雜變化指標的EMAP值相差較小，分別為2%和3%，截口指標的EMAP值分別為7%和14%，上述序列在復雜變化特性方面擬合效果相差較小，但保持實測序列截口統計特性的效果相差較大。單獨使用截口或復雜性指標進行評價，不能保證模擬序列年內分配方面的擬合效果。例如息縣站的2-2、2-10序列，分別在復雜性、截口指標中EMAP值較小，但其平、枯水年及多年平均情況下的相對誤差均較大。綜合上述分析，需要全面考慮模擬序列截口、復雜變化特性方面的擬合效果。

表4 模擬序列不同角度精度分析結果對比Tab.4 Comparison of accuracy analysis results of different angles of simulation sequence

3.4 模擬序列優選

基于截口-復雜性指標體系，采用灰關聯分析法對模擬序列進行優選，結果見表5。由表5可以看出漯河站、息縣站和蔣家集站最優模擬序列分別為1-5序列、2-1和3-7序列。

表5 模擬序列的關聯度Tab.5 Correlation degrees of simulated sequences

為進一步說明截口-復雜性指標體系用于模擬序列評價的合理性，對傳統指標、樣本熵最小指標[9]以及截口-復雜性指標體系優選出的模擬序列進行對比，截口、復雜變化特性方面的相對誤差和EMAP見表6和表7，多年平均以及豐、平、枯水年的擬合效果見圖2。

表6 3種評價指標優選序列截口指標的相對誤差及EMAP結果統計Tab.6 Statistics of RE and EMAP results of optimized sequence section indexes of three evaluation indexes %

表7 3種評價指標優選序列復雜性指標的相對誤差及EMAP結果統計Tab.7 Statistics of RE and EMAP results of complex change index of three evaluation index optimization series %

由表6和表7可以看出：3種指標優選的模擬序列截口指標的EMAP值差異較小，最大為2%；復雜變化指標的EMAP值差異較大，截口-復雜性指標體系和傳統指標優選序列的EMAP值最大相差8%，與樣本熵指標優選序列的EMAP值最大相差4%。由圖2可以看出，截口-復雜性指標體系優選序列的各水平年年內分配相對誤差總體小于傳統指標和樣本熵指標優選序列的結果，且平、枯水年相對誤差差異最為明顯。

圖2 3種評價指標優選序列多年平均及豐/平/枯水年相對誤差Fig.2 The relative errors of three evaluation index optimization series in annual average,wet,normal and dry years

綜合上述分析，3種指標評價方法優選序列的精度為：截口-復雜性指標>樣本熵指標>傳統指標。這說明本文提出的截口-復雜性指標體系是切實可行的，有利于全面評價徑流模擬序列，實現不同隨機模擬序列之間的比較與優選。

3.5 討 論

截口、復雜性及水平年內分配方面的對比結果表明，3種指標優選的序列均較好地保持了實測序列的截口特性，但傳統指標、樣本熵最小指標體系不能全面評價模擬序列，有可能遺漏徑流序列重要的細節信息。以下從年內分配異常情況和時程分配(過程)的差異性兩方面進行探討。

根據3個水文站歷年逐月實測徑流資料統計結果，漯河站、息縣站、蔣家集站多年平均汛期最大4個月(6—9月)徑流量占全年徑流量分別為67%、64%、61%。現將連續最大4個月徑流出現在非汛期視為異常情況，統計傳統截口指標優選序列、截口-復雜性指標體系優選序列年內連續最大4個月徑流在非汛期的百分率，并與實測資料進行對比(表8)，結果發現：傳統指標優選序列的百分率在3種序列類型中最大，且與截口-復雜性指標優選序列的情況相差較大。這說明只保證模擬系列均值、方差等基本統計特征較為一致情況下，模擬序列會出現較多的年內分配異常情況，而截口-復雜性指標體系有助于排除年內分配異常的情況。

表8 3種徑流序列連續最大4個月徑流在非汛期的百分率Tab.8 Percentage of the maximum 4-month continuous runoff of three runoff series in non-flood season

表9統計了漯河站3種序列的年內分配指標，可以看出：截口-復雜性指標體系優選序列的標準差數值整體大于傳統指標優選序列，說明本文指標體系增加的復雜變化指標并沒有影響模擬序列各子序列年內分配特征的多樣性；兩種指標優選序列與實測序列的標準差差異較小，說明傳統指標體系和截口-復雜性指標體系優選的長序列模擬序列中，各子序列的多樣性受多年平均計算條件下各指標約束的影響較小。

表9 3種徑流序列年內分配指標分析Tab.9 Analysis of annual distribution index of three runoff series

模擬序列截口、復雜變化特性的評價以及年內分配異常情況的檢驗表明，引入表征復雜變化特性的指標(Q4、Cd、Ct、H)，可以在保證模擬序列多樣性的同時，減少年內分配的異常情況，能夠描述傳統評價指標所忽略的細節特征，較好地體現原始序列的形態復雜變化特征。因此，截口-復雜性指標體系更準確、真實地刻畫了月徑流模擬序列的年內復雜變化特征，優選出的模擬序列更能保持實測序列的時間變化特征。

4 結 論

本文從徑流序列非線性及復雜特性出發，提出了月徑流隨機模擬序列評價指標體系，并將其用于3個水文站季節性自回歸模型的模擬序列效果評價中，得到以下結論。

只考慮截口統計參數的評價方法，忽視了徑流序列的復雜變化特性，可能遺漏重要的細節信息，從而影響模擬序列檢驗的精度。截口-復雜性指標體系中表征復雜變化特性的指標(Q4、Cd、Ct、H)能簡單有效地區分模擬序列保持實測資料復雜變化特性的能力，優選序列在年內分配方面的EMAP值最小，總體擬合效果最好。

對傳統指標、樣本熵最小指標以及截口-復雜性指標優選的序列進行了分析，結果表明截口-復雜性指標體系優選的序列能較好地反映實測序列的變化過程和復雜特性，綜合截口和復雜變化特征方面的分析，月徑流隨機模擬序列優選方法的優劣順序為截口-復雜性指標>樣本熵指標>傳統指標。

經過截口和復雜變化特征方面的模擬序列評價、年內分配異常情況分析以及時程分配(過程)差異性分析的多重驗證，截口-復雜性指標體系在保證模擬序列年內分配多樣性的同時可以減少年內分配異常的情況，有助于水文工作者掌握不同模擬序列之間的差異，更加全面地評價徑流模擬序列，從而提供更加準確的水文模擬序列。

為進一步提高徑流序列模擬效果評價的準確性，未來將考慮采用主成分分析法找出相互獨立的主要影響指標，以此構成更具代表性的指標體系，并由此對隨機模擬序列進行檢驗，通過與其他指標體系檢驗結果的對比分析其代表性。