拒馬河上游年徑流序列特性分析及其隨機模擬

王 玲

（河北省保定水文水資源勘測局，河北 保定071000）

1 基本情況

拒馬河是大清河的干流，主要為泉水，從淶源縣城發(fā)源，流經(jīng)易縣紫荊關(guān)，淶水野三坡，北京房山十渡，在房山張坊分為南拒馬河和北拒馬河。

北拒馬河流經(jīng)北京南尚樂鄉(xiāng)，于二合莊村東出市境，入河北涿州市境內(nèi)，至東茨村以下稱白溝河，在白溝村與南拒馬河匯合入大清河，匯入白洋淀［1］。

長期以來拒馬河上游未建大型水庫，河畔兩側(cè)多山巒，遇暴雨天氣，極易發(fā)生洪水，且來勢兇猛。紫荊關(guān)水文站位于拒馬河上游，對于拒馬河的防洪有著重要作用。2012年7月21～22日拒馬河發(fā)生了1963年以來最大的一場暴雨洪水，紫荊關(guān)水文站洪峰流量達2580m3/s，因此對紫荊關(guān)站徑流系列的隨機模擬具有重要意義。

紫荊關(guān)水文站位于河北省易縣境內(nèi)，流域分布在淶源縣、易縣境內(nèi)，地勢西高東低，山巒起伏為太行山脈，一般山峰海拔在2100m以下，紫荊關(guān)水文站流域面積1760km2，主河道長81.5km，河道縱坡5.5%，流域平均寬度25.4km，最大洪峰流量4490m3/s（1963年）［2］。

2 資料審查

選用紫荊關(guān)站徑流序列資料 （1965～2013年）均系水文部門實測資料。多年來，位于拒馬河上游的紫荊關(guān)站集水面積和氣候條件并沒有發(fā)生大的變化，也沒有大的水利工程建設(shè)，從而能夠滿足資料一致性的要求。由于紫荊關(guān)水文站40多年的年徑流序列有很好的穩(wěn)定性，而且表現(xiàn)為豐枯水年組的周期性變化，可以表明所選徑流資料有很好的代表性。

3 紫荊關(guān)站年徑流序列自相關(guān)性分析及模型識別

根據(jù)模型識別的基本原理，徑流序列究竟符合哪一類隨機模擬模型，取決于該徑流序列的自相關(guān)性和其偏態(tài)特性，為此首先分析紫荊關(guān)站年徑流序列的自相關(guān)性［3］。

年徑流序列Xt的自相關(guān)系數(shù)表示為：

式中 rk為年徑流序列Xt的自相關(guān)系數(shù)；n為序列長度，當(dāng)n＞50時，可以取m＜n/4，常取m在n/10左右；當(dāng)n＜50時，取m在n/4左右的數(shù)值，有的取m＜n-10；k為年徑流序列的滯時；為年徑流序列的均值。

本文所選紫荊關(guān)站年徑流序列年數(shù)為49a，取m=12，計算自相關(guān)系數(shù)并繪制自相關(guān)圖，并加繪了獨立序列自相關(guān)置信水平為95%的容許限，如表1和圖1。

表1 紫荊關(guān)站49a徑流序列自相關(guān)系數(shù)及其容許限

圖1 紫荊關(guān)站年徑流自相關(guān)圖

表1和圖1表明：紫荊關(guān)站年徑流量序列的一階和二階自相關(guān)系數(shù)均顯著異于獨立序列，因此為一組相依序列。對正相依序列，可以選用不少類型的模型來描述其統(tǒng)計特征。考慮以下幾點，選用自回歸模型AR（p）來模擬紫荊關(guān)站的徑流序列。

（1）AR（p）模型表征徑流序列的統(tǒng)計特性有一定的物理基礎(chǔ)；

（2）AR（p）模型參數(shù)的估計可以用簡單的矩法，而且精度高；

（3）AR（p）模型形式簡單，數(shù)學(xué)處理方法簡便，為大家熟悉。

4 紫荊關(guān)站年徑流序列偏態(tài)特性分析

選定AR（p）模型后，如何確定p，需分析其序列的偏態(tài)特性，主要看是否在k＞p后，φp，p偏態(tài)系數(shù)在95%的容許限范圍內(nèi)，當(dāng)假設(shè)成立，可以認為能夠用AR（k）模擬該徑流序列［4］。

對于自回歸AR（p）模型有

式中 φp，1為p階模型的第一系數(shù)；φp，p為p階模型的第p個系數(shù)。

式（6）中，偏相關(guān)系數(shù)φp，p是反映消除（p-1）階自相關(guān)系數(shù)影響后所剩余的自相關(guān)程度，計算方法上，偏相關(guān)函數(shù)和自相關(guān)函數(shù)有一定聯(lián)系，利用自相關(guān)函數(shù)可以計算偏相關(guān)函數(shù)。可令φk，j表示自回歸AR（k）模型中的第j個系數(shù)，這樣φk，k為最后一個系數(shù)。則：

可求的φk，j（j=1，2，…，k），對于自回歸AR（1）模型

對于高階自回歸模型，可應(yīng)用以下述遞推算法：

計算結(jié)果如表2，偏相關(guān)函數(shù)圖如圖2。

表2 紫荊關(guān)站49a徑流序列偏相關(guān)系數(shù)及其容許值

續(xù)表2

圖2 紫荊關(guān)站年徑流偏相關(guān)圖

圖2表明：當(dāng)k≥2時，φk，k落于容許限內(nèi)，故可推斷p=1。換言之，據(jù)偏相關(guān)系數(shù)的統(tǒng)計分析，AR（1）模型可以用來描述年徑流量序列的統(tǒng)計變化。但是，在序列長度n相當(dāng)大且序列為正態(tài)分布的情況下才是完全正確的，就紫荊關(guān)年徑流量序列而言，其長度僅有49a，而其二階自相關(guān)系數(shù)較大，這暗示該序列可能為AR（2）序列。有以上分析，可以認為紫荊關(guān)年徑流量序列可能是AR（1）序列，也可能是AR（2）序列。下面對這兩種可能的模型均作參數(shù)估計，以便進一步檢驗。

5 參數(shù)估計及年徑流序列模擬

5.1 參數(shù)估計

通過對參數(shù)的平穩(wěn)性分析，最終求得AR（1）模型和AR（2）模型。

AR（1）模型為：

式中 φt為服從均值0、方差1的偏態(tài)系數(shù)；Csφ的P-Ⅲ型分布，用公式（11），（12）可求得。

AR（2）模型為：

5.2 模型形式的進一步判斷

估計出參數(shù)后，便可利用AIC準(zhǔn)則進一步識別紫荊關(guān)站年徑流序列是AR（1）還是AR（2）模型。由式（11）計算AIC值，從表3可以看出，由AIC準(zhǔn)則來判斷，AR（1）模型優(yōu)于AR（2）模型。

表3 紫荊關(guān)站49a徑流序列AIC值

5.3 模型檢驗

主要檢驗εt是否獨立，對AR （1）模型，令εt=0.735φt，則有：

根據(jù)實測樣本序列可計算得εt，由此計算出r1（ε），r2（ε），…，r12（ε），最后計算統(tǒng)計量Q=8.91。根據(jù)自由度和顯著水平，得Q<χ2α，故為獨立的假設(shè)可以接受。

φt為自回歸模型的隨機變量部分，用式（11）、式（12）計算：

式中 Csφ為φ分布的偏態(tài)系數(shù)；Csx為X分布的偏態(tài)系數(shù)；εt為隨機數(shù)。

模擬隨機數(shù)εt時，本文采用變換法，即對均勻隨機數(shù)做下列變換：

則ξ1，ξ2為相互獨立的標(biāo)準(zhǔn)正態(tài)分布N（0，1）變量，依次模擬出ξt。

6 模型檢驗和實用性分析

重點分析AR（1）能否保持實測序列的主要統(tǒng)計特征［5］。分別根據(jù)AR（1）模型的模擬徑流序列，然后按長序列和短序列法計算各種參數(shù)，并和實測序列的相應(yīng)參數(shù)做對比，結(jié)果如表4。

表4 模擬序列的實用性檢驗

由表4可以看出，長序列的模擬結(jié)果優(yōu)于短序列的結(jié)果，長序列的AR（1）的基本參數(shù)的誤差在10%以內(nèi)，均能很好地保持，短序列模擬的各參數(shù)均偏小，這種偏小是由于計算Cv，Cs，r等公式在樣本容量較小時皆為負偏而造成。因此使用一階自回歸模擬AR（1）能夠較好地模擬拒馬河上游的年徑流序列，基本滿足實際工程應(yīng)用需要。

7 結(jié)語

通過對拒馬河水系上游紫荊關(guān)站年徑流變化趨勢分析和隨機模擬，可以清晰掌握拒馬河的年徑流的變化規(guī)律。通過進一步的水文模擬，可以更清楚地看出拒馬河紫荊關(guān)站年徑流量的變化范圍，從而能夠在一定程度上服務(wù)于紫荊關(guān)站的防洪和水利工程建設(shè)。

［1］郭慶宏，朱虹.河北省拒馬河山區(qū)降水特性分析［J］.河北工程技術(shù)高等專科學(xué)校學(xué)報，2012（2）：1-4.

［2］李紹飛，余萍，等.紫荊關(guān)流域洪水徑流過程變化及影響因素分析［J］.武漢大學(xué)學(xué)報（工學(xué)版），2012，45（2）：167-170.

［3］丁晶，鄧育仁.隨機水文學(xué)［M］.成都：成都科技大學(xué)出版社，1998.

［4］Bras R A.Random Functions and Hydrology ［M］.America：Addison Wesley Publication Co.，1985.

［5］金光炎.水文統(tǒng)計原理與方法［M］.北京：中國工業(yè)出版社，1954.