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

王 玲

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

1 基本情況

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

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

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

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

2 資料審查

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)于自回歸AR（p）模型有

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

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

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

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

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

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

5.1 參數(shù)估計(jì)

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

AR（1）模型為：

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

AR（2）模型為：

5.2 模型形式的進(jìn)一步判斷

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

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

5.3 模型檢驗(yàn)

主要檢驗(yàn)εt是否獨(dú)立，對(duì)AR （1）模型，令εt=0.735φt，則有：

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

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

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

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

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

6 模型檢驗(yàn)和實(shí)用性分析

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

表4 模擬序列的實(shí)用性檢驗(yàn)

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

7 結(jié)語(yǔ)

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

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

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

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

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

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