基于ARMA模型的水文序列相依變異分級(jí)方法及驗(yàn)證

謝 平，霍競(jìng)?cè)海Ｑ喾?，吳林倩，李雅晴，牛靜怡

（1.武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430072；2.中國(guó)科學(xué)院地理科學(xué)與資源研究所陸地水循環(huán)與地表過(guò)程重點(diǎn)實(shí)驗(yàn)室，北京 100101）

1 研究背景

水文時(shí)間序列分析是深入了解水文現(xiàn)象特點(diǎn)與性質(zhì)的重要途徑，以及探究水文過(guò)程復(fù)雜變化規(guī)律的關(guān)鍵方法［1-2］。長(zhǎng)久以來(lái)，對(duì)水文時(shí)間序列的認(rèn)識(shí)與研究都是基于一致性假設(shè)（即物理成因不變）的基本條件開(kāi)展的［1］。然而，受自然因素與人為因素等的影響，許多流域和地區(qū)的水文水資源狀況發(fā)生了很大改變，導(dǎo)致水文時(shí)間序列不再全部滿足一致性假設(shè)［3-4］，即發(fā)生了水文變異［5］，表現(xiàn)為統(tǒng)計(jì)規(guī)律隨時(shí)間發(fā)生明顯改變，且任意時(shí)刻數(shù)值均與前期時(shí)間內(nèi)的歷史值有關(guān)［6］。

依據(jù)水文變異的具體表現(xiàn)形式與水文序列參數(shù)的變化特征，又可分為跳躍、趨勢(shì)、周期以及相依變異等［7］。其中，相依變異是水文變異的一種常見(jiàn)形式，表現(xiàn)為水文序列后一數(shù)據(jù)與緊鄰的前一（幾）數(shù)據(jù)具有成因、數(shù)值等方面的繼承性［2］。當(dāng)前，對(duì)水文序列中相依變異成分的研究和分析主要是基于數(shù)理統(tǒng)計(jì)方法進(jìn)行［8-9］，如Hurst系數(shù)法［10-11］、極差和輪次分析法［12］、自相關(guān)分析法和譜分析法［13］等。然而，這些方法主要針對(duì)相依變異形式的識(shí)別和檢驗(yàn)，但對(duì)相依變異程度進(jìn)行量化的研究則很少涉及。在實(shí)際的水文統(tǒng)計(jì)分析與計(jì)算中，不僅要明確水文序列是否發(fā)生了相依變異，還要判定相依變異的強(qiáng)弱程度，如果未對(duì)相依變異的顯著性進(jìn)行定量評(píng)估，可能會(huì)對(duì)水資源安全評(píng)估以及水利工程建設(shè)等工作造成較大影響。相關(guān)系數(shù)可以描述序列間相關(guān)程度，其絕對(duì)值越大，則序列間的相關(guān)性越強(qiáng)，故可將相關(guān)系數(shù)作為度量水文變異程度的指標(biāo)［14］。

常規(guī)的時(shí)間序列模型包括自回歸模型、滑動(dòng)平均模型和自回歸滑動(dòng)平均模型。因其具有時(shí)間相依的直觀形式和較強(qiáng)的適用性，故水文序列的相依性一般用它們來(lái)進(jìn)行描述［15］。其中，自回歸滑動(dòng)平均模型ARMA（Auto-Regressive Moving Average）是自回歸模型AR（Auto-Regressive）和滑動(dòng)平均模型MA（Moving Average）的混合，且較AR和MA能更好地反映水文變量在時(shí)序變化上的統(tǒng)計(jì)特性，即具有更大的彈性，在水文學(xué)中應(yīng)用更為廣泛［7］。國(guó)內(nèi)外許多學(xué)者用其來(lái)描述水文現(xiàn)象在時(shí)間上的相依性。翟顥瑾等［16］使用ARMA模型模擬水質(zhì)污染的相依隨機(jī)過(guò)程，實(shí)現(xiàn)對(duì)長(zhǎng)江未來(lái)水質(zhì)情況的分析預(yù)測(cè)；栗現(xiàn)文和楊佳等［17-18］采用ARMA模型分析地下水位埋深時(shí)間序列中去除確定性組分之后的相依隨機(jī)和純隨機(jī)組分，從而達(dá)到對(duì)未來(lái)短期地下水位動(dòng)態(tài)進(jìn)行準(zhǔn)確預(yù)測(cè)的目的；Dwivedi等［19］運(yùn)用ARMA等時(shí)間序列模型對(duì)降雨和徑流的隨機(jī)相依水文過(guò)程建模，從而實(shí)現(xiàn)對(duì)降雨的預(yù)測(cè)和徑流量的估計(jì)；Atan 等［20］利用ARMA 模型的相依性特征模擬洪水時(shí)間序列并對(duì)流域洪災(zāi)風(fēng)險(xiǎn)進(jìn)行了評(píng)估與分析等。然而，目前的研究主要是運(yùn)用ARMA模型表征水文序列的相依特性，對(duì)于水文序列中相依變異強(qiáng)弱的量化研究則鮮有提及。趙羽西等［21］的研究雖涉及到對(duì)序列相依變異程度的分級(jí)量化，卻也僅適用于AR 模型，至于其它模型尚需進(jìn)一步補(bǔ)充和完善。此外，在運(yùn)用AIC、BIC 等定階準(zhǔn)則對(duì)含有AR?MA相依成分的水文序列進(jìn)行定階時(shí)，只能確定模型的自回歸階數(shù)p與滑動(dòng)平均階數(shù)q的和（即p+q的值），而在AIC和BIC函數(shù)最小值不唯一或階數(shù)p和q有多種組合方式時(shí)，并不能對(duì)其做出明確判斷。

為此，本文以適用性更為廣泛的ARMA模型為例，在已有研究［21］的基礎(chǔ)上，以原序列與其相依成分間的相關(guān)系數(shù)為衡量標(biāo)準(zhǔn)，提出了描述水文相依變異強(qiáng)弱的一種方法。通過(guò)相關(guān)系數(shù)表達(dá)式的推求說(shuō)明選用該指標(biāo)進(jìn)行變異分級(jí)的基本原理，并設(shè)計(jì)統(tǒng)計(jì)試驗(yàn)說(shuō)明基本原理的合理性；以ARMA（1，1）和ARMA（1，2）模型為例，指出相關(guān)系數(shù)指標(biāo)與自相關(guān)系數(shù)的聯(lián)系；最后通過(guò)模擬時(shí)間序列和實(shí)測(cè)水文序列對(duì)所提方法進(jìn)行驗(yàn)證，并結(jié)合水文序列的相依變異程度和相依成分與原序列的擬合效率系數(shù)為具有ARMA（p，q）相依成分的水文序列提供一種定階的新思路。

2 分級(jí)原理與方法

2.1 相依變異分級(jí)原理一般情況下，在扣除確定性成分的影響之后，剩余水文序列xt中的相依成分可用自回歸滑動(dòng)平均模型ARMA（p，q）表示為：

式中：u為序列xt的均值；p為模型的自回歸階數(shù)；q為模型的滑動(dòng)平均階數(shù)；參數(shù)φ1，φ2，…，φp是待定系數(shù)中的自回歸系數(shù)；參數(shù)θ1，θ2，…，θq是待定系數(shù)中的滑動(dòng)平均系數(shù)；εt，εt-1，εt-2，…，εt-q為獨(dú)立隨機(jī)變量，其標(biāo)準(zhǔn)差為σε、均值為0，且與xt-1，xt-2，…，xt-q獨(dú)立無(wú)關(guān)。

式（1）中的φ1(xt-1-u)+φ2(xt-2-u)+…+φp(xt-p-u)-θ1εt-1-θ2εt-2-…-θqεt-q為相依成分，用ηt表示，故xt能夠表示為線性疊加的形式：

式中：ut=u+εt，其均值與標(biāo)準(zhǔn)差分別為u、σε，為水文序列中的純隨機(jī)成分，在我國(guó)水文統(tǒng)計(jì)分析與計(jì)算中通常認(rèn)為其服從P-Ⅲ型分布［7］。

又因?yàn)棣莟與ut相互獨(dú)立，而xt的期望E(xt)=u、ut的期望E(ut)=u，故可得相依成分序列的期望值E(ηt)=0。

xt與ηt間的相關(guān)系數(shù)表達(dá)式為：

式中：

因此，由式（4）可得：

式中ση2與σx2分別為ηt和xt的方差，它們之間具有以下關(guān)系：

式中σu2為純隨機(jī)序列的方差，根據(jù)式（6）和式（7）可得：

將式（1）進(jìn)行移項(xiàng)并整理得：

將式（9）兩邊同乘以xt-u并取數(shù)學(xué)期望得：

上式兩端除以σx2后考慮到自相關(guān)系數(shù)的定義：

可得：

式（12）兩端再同乘以σx2可得：

將式（9）兩邊同乘以εt并取數(shù)學(xué)期望得：

將式（9）兩邊同乘以εt-1并取數(shù)學(xué)期望得：

同理，依此類(lèi)推，直至將式（9）兩邊同乘以εt-q再取數(shù)學(xué)期望得：

為此，式（13）可以化簡(jiǎn)為：

將上式進(jìn)行整理得：

因?yàn)棣姚?=σu2，結(jié)合式（8）和（18）可得：

對(duì)式（19）進(jìn)行初步驗(yàn)證：

當(dāng)p=0時(shí)，即為滑動(dòng)平均模型MA（q），式（19）可化簡(jiǎn)為：

當(dāng)q=0時(shí)，即為自回歸模型AR（p），式（19）可化簡(jiǎn)為：

式（20）和（21）與現(xiàn)有結(jié)論相符合［22］，故初步認(rèn)為式（19）合理可靠。

ARMA 模型中的滑動(dòng)平均系數(shù)θ1，θ2，…，θq和自回歸系數(shù)φ1，φ2，…，φp的參數(shù)解與樣本自相關(guān)系數(shù)存在函數(shù)關(guān)系，而樣本自相關(guān)系數(shù)能夠描述水文序列的相依變異程度［7］，故相關(guān)系數(shù)r能夠表示樣本序列的相依關(guān)系并可描述其相依程度。因此，可以用相關(guān)系數(shù)對(duì)水文序列的相依程度進(jìn)行分級(jí)。

2.2 相依變異分級(jí)方法基于以上推導(dǎo)證明的分級(jí)原理，借鑒水文變異分級(jí)方法［21，23-25］，提出水文時(shí)間序列的相依變異分級(jí)方法。其中，選取rα、rβ、0.6和0.8作為相依變異分級(jí)界限，相應(yīng)地把相關(guān)系數(shù)r分為5段區(qū)間，對(duì)應(yīng)描述不同強(qiáng)弱的相依變異程度，具體分級(jí)標(biāo)準(zhǔn)見(jiàn)表1。rα和rβ分別是在已知模型階數(shù)和序列長(zhǎng)度的情況下，顯著性水平為α和β（α>β）下的分級(jí)閾值，通常取α=0.05，β=0.01。

表1 相依變異程度分級(jí)

以ARMA（p，q）為例，利用所提方法對(duì)水文序列的相依變異程度量化分級(jí)的詳細(xì)步驟如下：

（1）去除水文序列中的確定性成分，得到含有獨(dú)立隨機(jī)或相依成分的剩余序列xt。

（2）判斷xt中是否存在相依成分。畫(huà)出自相關(guān)與偏相關(guān)系數(shù)圖，如果自相關(guān)系數(shù)均在容許限以內(nèi)，認(rèn)為剩余序列為純隨機(jī)成分；否則，認(rèn)為剩余序列含有相依成分，且當(dāng)自相關(guān)與偏相關(guān)系數(shù)均拖尾時(shí)，認(rèn)為相依成分符合ARMA（p，q）模型。

（3）確定模型的階數(shù)。不同于AR（p）與MA（q），據(jù)自相關(guān)系數(shù)圖和偏相關(guān)系數(shù)圖并不能判斷出ARMA（p，q）模型的階數(shù)，故選用兩種定階準(zhǔn)則（AIC和BIC準(zhǔn)則）來(lái)判定其階數(shù)p、q。

（4）估計(jì)模型參數(shù)。用矩估計(jì)法［26］對(duì)AR?MA（p，q）的參數(shù)φi(i=1，2，...，p)和θi(i=1，2，...，q)進(jìn)行估計(jì)，由此得到相依成分序列ηt。

（5）計(jì)算相關(guān)系數(shù)r。利用式（3）求得序列xt與ηt間的r值。

（6）判定變異等級(jí)。根據(jù)得出的相關(guān)系數(shù)r值對(duì)照表1找到對(duì)應(yīng)的等級(jí)區(qū)間，由此判定出水文序列相依變異程度的強(qiáng)弱。

3 分級(jí)方法驗(yàn)證

3.1 分級(jí)原理合理性驗(yàn)證為驗(yàn)證水文相依變異分級(jí)方法對(duì)ARMA（p，q）的適用性，而實(shí)際中基于時(shí)間序列模型的水文序列預(yù)測(cè)與延展工作通常以低階為主，故以較低階數(shù)（p+q）的ARMA（1，1）和ARMA（1，2）模型為例，分別設(shè)計(jì)統(tǒng)計(jì)試驗(yàn)，以驗(yàn)證分級(jí)原理中式（19）推導(dǎo)過(guò)程的合理性。

3.1.1ARMA（1，1）模型ARMA（1，1）模型可以表示為xt=φ1(xt-1-u)-θ1εt-1+ut，據(jù)式（19）可得：

用Yule-Walker方程估計(jì)自回歸系數(shù)φi(i=1，2，…，p)［27］：

由Yule-Walker方程可得ρ1=φ1，將其代入式（22）得到：

為滿足平穩(wěn)性和可逆性，各系數(shù)需符合：

統(tǒng)計(jì)試驗(yàn)詳細(xì)步驟如下：

（1）利用P-Ⅲ型分布生成長(zhǎng)為n=1000、均值u=100、變差系數(shù)Cvu=0.2、偏態(tài)系數(shù)Csu=0.4 的序列ut，然后隨機(jī)生成20組自回歸系數(shù)φ1值和滑動(dòng)平均系數(shù)θ1值，代入式（1）得到20個(gè)ARMA（1，1）序列xi（i=1，2，3，…，20）；

（2）相依成分ηi由所得ARMA（1，1）序列xi扣除ut求??；

（3）每組系數(shù)值內(nèi)試驗(yàn)1000 次，ηij為第i組第j次統(tǒng)計(jì)試驗(yàn)得到的相依成分序列（i=1，2，…，20；j=1，2，…，1000）；

（4）依據(jù)式（3）計(jì)算得到每組原始序列與其相依成分序列的相關(guān)系數(shù)rij（i=1，2，···，20；j=1，2，…，1000），再求出每組的平均相關(guān)系數(shù)

表2 ARMA（1，1）模型中不同參數(shù)下的ra、ri 及相對(duì)誤差（δ）

3.1.2ARMA（1，2）模型ARMA（1，2）模型可以表示為xt=φ1(xt-1-u)-θ1εt-1-θ2εt-2+ut，將p=1，q=2 帶入式（19）得：

因?yàn)橛蒠ule-Walker方程得到ρ1=φ1，代入式（26）得：

為滿足平穩(wěn)性和可逆性，各系數(shù)需符合如下條件：

ARMA（1，2）模型的統(tǒng)計(jì)試驗(yàn)步驟參照ARMA（1，1）模型，并保持參數(shù)不變。在符合式（28）要求的情況下，隨機(jī)生成20個(gè)不同φ1、θ1、θ2的系數(shù)值組合，依據(jù)模型表達(dá)式得到ARMA（1，2）序列，對(duì)比統(tǒng)計(jì)試驗(yàn)得到的平均相關(guān)系數(shù)ri與由公式求出的相關(guān)系數(shù)ra，結(jié)果見(jiàn)表3?？梢钥闯?，20組統(tǒng)計(jì)試驗(yàn)中只有兩組結(jié)果的δ值超出5%，但不足6%，且參數(shù)個(gè)數(shù)越少，相對(duì)誤差δ越小，由此驗(yàn)證了式（19）的合理性，說(shuō)明本文所提分級(jí)原理對(duì)ARMA（1，2）模型也具有適用性。

表3 ARMA（1，2）模型中不同參數(shù)下的ra、ri 及相對(duì)誤差（δ）

此外，進(jìn)行了針對(duì)ARMA（2，1）、ARMA（1，3）、ARMA（2，2）和ARMA（3，1）模型的統(tǒng)計(jì)試驗(yàn)，且根據(jù)統(tǒng)計(jì)試驗(yàn)結(jié)果，本文所提分級(jí)原理與方法對(duì)于其它階數(shù)p+q等于4階以內(nèi)的ARMA模型也是合理的，限于篇幅，詳細(xì)分析結(jié)果不再展開(kāi)贅述。因此，本文所提的分級(jí)原理與方法對(duì)較低階數(shù)的AR?MA（p，q）模型具有普適性。

3.2 相關(guān)系數(shù)與自相關(guān)系數(shù)的關(guān)系式（19）給出了r與自相關(guān)系數(shù)的關(guān)系，其成為本文所提分級(jí)方法的理論基礎(chǔ)。為使該分級(jí)方法的合理性更具說(shuō)服力，現(xiàn)繼續(xù)探討相關(guān)系數(shù)與自相關(guān)系數(shù)在ARMA模型中的具體聯(lián)系，此處僅以ARMA（1，1）和ARMA（1，2）模型為例。

3.2.1ARMA（1，1）模型 利用統(tǒng)計(jì)試驗(yàn)法生成具有ARMA（1，1）相依成分的序列，隨機(jī)序列ut服從P-Ⅲ型分布，其長(zhǎng)度n=1000、均值u=100、變差系數(shù)Cvu=0.2、偏態(tài)系數(shù)Csu=0.4。為單獨(dú)研究?jī)煞N待定系數(shù)的影響，選取兩組φ1、θ1值。第一組滑動(dòng)平均系數(shù)θ1取固定值0.2，自回歸系數(shù)分別取φ1=-0.9、-0.5和-0.2 ； 第二組自回歸系數(shù)φ1取固定值0.2， 滑動(dòng)平均系數(shù)分別取θ1=0.2、0.5和0.9。圖1、圖2為兩組情況所得各階自相關(guān)系數(shù)。

圖1和圖2顯示， |ρ1|不僅隨著 |φ1|的增大而增大，也隨著 |θ1|增加而增加。而且，序列與相依成分的相關(guān)系數(shù)r也存在與 |ρ1|一致的變化趨勢(shì)，即隨著 |φ1|與 |θ1|的增大而增大。由此可知，對(duì)含有ARMA（1，1）相依成分的序列而言，相關(guān)系數(shù)可以刻畫(huà)其相依變異程度。

圖1 滑動(dòng)平均系數(shù)θ1 固定時(shí)的各階自相關(guān)系數(shù)

圖2 自回歸系數(shù)φ1固定時(shí)的各階自相關(guān)系數(shù)

3.2.2ARMA（1，2）模型 自相關(guān)系數(shù)可以描述水文序列的相依程度［7］，在ARMA（1，2）模型中，式（27）顯示r2與參數(shù)φ1、θ1、θ2有關(guān)，而滑動(dòng)平均系數(shù)和自回歸系數(shù)的參數(shù)解與樣本自相關(guān)系數(shù)存在函數(shù)關(guān)系，所以r與自相關(guān)系數(shù)具有函數(shù)關(guān)系。根據(jù)式（28）以及自相關(guān)系數(shù)、自回歸系數(shù)和滑動(dòng)平均系數(shù)的關(guān)系，選取多組符合條件的ρ1、ρ2的值，并由矩估計(jì)法［26］得到ARMA（1，2）模型的參數(shù)φ1、θ1、θ2值；純隨機(jī)成分的生成同3.2.1 節(jié)。ARMA（1，2）序列減去純隨機(jī)成分得到相依成分序列。求出每種ρ1和ρ2值下的相關(guān)系數(shù)，結(jié)果如表4所示。能夠看出：一階自回歸系數(shù)φ1的大小與ρ1的值無(wú)關(guān)［26，28］，且當(dāng)二階自相關(guān)系數(shù)ρ2固定不變時(shí)，相關(guān)系數(shù)r隨著一階自相關(guān)系數(shù)ρ1的增加而增加；反之，當(dāng)ρ1固定時(shí)，相關(guān)系數(shù)r隨著ρ2的增加而減小。說(shuō)明r與 |ρ1|呈正相關(guān)關(guān)系、與 |ρ2|呈負(fù)相關(guān)關(guān)系，故對(duì)于含有ARMA（1，2）相依成分的序列，其相依性強(qiáng)弱也能用相關(guān)系數(shù)進(jìn)行描述。

表4 不同ρ1、 ρ2 下序列的相關(guān)系數(shù)

3.3 模擬序列分析驗(yàn)證將所提分級(jí)方法應(yīng)用于具有相依成分的模擬序列進(jìn)行驗(yàn)證。首先從常用的顯著性水平中選取第一顯著性水平α=0.05、第二顯著性水平β=0.01；然后根據(jù)相依變異程度分級(jí)表（即表1）中各個(gè)變異等級(jí)對(duì)應(yīng)的相關(guān)系數(shù)r的取值范圍，并結(jié)合所提水文相依變異分級(jí)方法的理論基礎(chǔ)式（19），求取各相依變異程度所對(duì)應(yīng)的整體變量的取值區(qū)間；最后選取滿足整體變量H在取值區(qū)間內(nèi)的系數(shù)值組合，生成具有相依變異成分的模擬序列。下面仍以AR?MA（1，1）和ARMA（1，2）模型為例進(jìn)行驗(yàn)證。純隨機(jī)成分在模擬序列中的長(zhǎng)度n=100、均值u=100、變差系數(shù)Cvu=0.2、偏態(tài)系數(shù)Csu=0.4，且服從P-Ⅲ型分布。

3.3.1ARMA（1，1）模型 在ARMA（1，1）模型中，它的自回歸部分和滑動(dòng)平均部分的階數(shù)分別為p=1、q=1，根據(jù)式（22）—（24）可知此時(shí)整體變量H可表示為。各相依變異程度對(duì)應(yīng)的整體變量H的取值區(qū)間見(jiàn)表5。

表5 ARMA（1，1）模型各相依變異程度區(qū)間所對(duì)應(yīng)的整體變量H的取值范圍

在上述條件下，分別選取一階自回歸與一階滑動(dòng)平均φ1=0.12、θ1=0.10，φ1=0.20、θ1=0.16，φ1=0.25、θ1=0.24，φ1=0.50、θ1=0.70 和φ1=0.80、θ1=0.95 共5 組系數(shù)值，由這5 組系數(shù)值求得的整體變量H分別對(duì)應(yīng)相依變異程度的5個(gè)等級(jí)。模擬生成這5種情形下的序列，繪于圖3，并計(jì)算出所生成序列與其相依成分之間的相關(guān)系數(shù)。

圖3為相依無(wú)、弱、中、強(qiáng)、巨變異下模擬序列與其相依成分示意圖。從圖3可以直觀看出，各相依成分與原序列的擬合程度依次加強(qiáng)，表明模擬序列的相依變異程度依次加強(qiáng)；求得各圖的模擬序列與相依成分間的相關(guān)系數(shù)值分別為0.1556、0.2530、0.3370、0.7047、0.9004，滿足表5所列出的各級(jí)相依變異程度對(duì)應(yīng)的相關(guān)系數(shù)的取值區(qū)間，也從定量角度很好地驗(yàn)證了所提方法的合理性。

圖3 ARMA（1，1）模型不同H值下相依變異分級(jí)序列

3.3.2ARMA（1，2）模型ARMA（1，2）模型的階數(shù)為p=1、q=2，由式（26）、（27）可知整體變量H可表示為，各相依變異程度所對(duì)應(yīng)的整體變量H的取值范圍見(jiàn)表6。

表6 ARMA（1，2）模型各相依變異程度區(qū)間所對(duì)應(yīng)的整體變量H的取值范圍

此種條件下，分別在相依變異的5 個(gè)等級(jí)對(duì)應(yīng)的整體變量H區(qū)間內(nèi)各取一組滿足條件的φ1、θ1、θ2組合值，同時(shí)一階自回歸系數(shù)φ1、一階滑動(dòng)平均系數(shù)θ1以及二階滑動(dòng)平均系數(shù)θ2還需滿足式（28）的限制條件。選取的5 組系數(shù)值分別為φ1=0.02、θ1=0.07、θ2=0.08，φ1=0.15、θ1=0.18、θ2=0.20，φ1=0.30、θ1=0.28、θ2=0.30，φ1=0.50、θ1=-0.50、θ2=-0.40，φ1=0.80、θ1=-0.80、θ2=-0.85。所生成序列與其相依成分的示意圖如圖4所示，圖4依次為相依變異從弱到強(qiáng)的5 個(gè)等級(jí)，所對(duì)應(yīng)的相關(guān)系數(shù)r值分別為0.1077、0.2975、0.4703、0.6842、0.9207，均滿足表6中各變異等級(jí)對(duì)應(yīng)的r值區(qū)間。圖4直觀地顯示了模擬序列相依變異的程度依次增強(qiáng)，而相關(guān)系數(shù)r的值也依次增加，故所提分級(jí)方法的合理性得到了驗(yàn)證。

圖4 ARMA（1，2）模型不同H值下相依變異分級(jí)序列

4 實(shí)測(cè)序列分析

為進(jìn)一步說(shuō)明所提分級(jí)方法的合理性與適用性，選取實(shí)測(cè)水文序列對(duì)其進(jìn)行驗(yàn)證，包括黃河花園口站的年徑流序列（1919—2000年）、長(zhǎng)江漢口（武漢關(guān)）站的月徑流序列（2010—2015年）及長(zhǎng)江沙市（二郎磯）站的月徑流序列（1991—1998年）。按照相依變異分級(jí)方法的步驟順序，要先扣除水文序列中的確定性成分，得到含有獨(dú)立隨機(jī)或相依成分的剩余序列，可以采用綜合加權(quán)法［1］、滑動(dòng)相關(guān)系數(shù)法［29］、Brown-Forsythe 法［30］等檢測(cè)序列的跳躍點(diǎn)；采用分段趨勢(shì)相關(guān)系數(shù)識(shí)別法［31］、過(guò)程線法、Kendall檢驗(yàn)法［32］等鑒別趨勢(shì)成分；采用滑動(dòng)周期相關(guān)系數(shù)識(shí)別法［33］、周期圖法、最大熵譜分析法［34-35］等檢驗(yàn)序列的周期成分。

顯著性水平分別選取α=0.05、β=0.01，對(duì)以上序列進(jìn)行變異診斷分析［36］，其中結(jié)果如表7所示。結(jié)果表明，花園口站的年徑流序列在1932年發(fā)生了明顯的跳躍變異，這主要是氣候變異直接導(dǎo)致的［37］；漢口站的月徑流序列和沙市站的月徑流序列都存在12 a 的周期變異，這種周期波動(dòng)現(xiàn)象主要和區(qū)域的氣候條件及水文地理特征密切相關(guān)［24］。

表7 實(shí)測(cè)徑流序列變異診斷結(jié)果

圖5為各站點(diǎn)實(shí)測(cè)徑流去除確定性成分之后的自相關(guān)系數(shù)和偏相關(guān)系數(shù)圖，藍(lán)線為容許度為95%時(shí)的界限。由圖可知，花園口站的自相關(guān)系數(shù)均在容許限以內(nèi)，說(shuō)明花園口站年徑流序列不含有相依成分；而漢口站和沙市站序列的自相關(guān)系數(shù)和偏相關(guān)系數(shù)均呈現(xiàn)出拖尾的特性，表明漢口站和沙市站的月徑流序列都存在相依成分，且為ARMA（p，q）相依成分。故采用ARMA（p，q）模型對(duì)漢口站和沙市站的月徑流的剩余序列進(jìn)行建模，ARMA（p，q）中的階數(shù)p，q運(yùn)用AIC與BIC準(zhǔn)則進(jìn)行判定：

圖5 實(shí)測(cè)水文序列的自相關(guān)系數(shù)及偏相關(guān)系數(shù)

式中：n與σε2分別為序列長(zhǎng)、殘差方差。AIC準(zhǔn)則和BIC準(zhǔn)則類(lèi)似，都是確定在有最小函數(shù)值時(shí)所對(duì)應(yīng)的p和q為最佳階數(shù)。

可得在p+q等于3和2時(shí)，漢口站和沙市站有最小的AIC與BIC值，由此可確定沙市站月徑流的剩余序列符合ARMA（1，1）模型，因在ARMA（p，q）模型中q的取值一般小于p［28］，故初步判斷漢口站月徑流的剩余序列符合ARMA（2，1）模型，還需要根據(jù)相依成分序列與其原始水文序列的擬合效果進(jìn)一步進(jìn)行驗(yàn)證。采用ARMA（2，1）和ARMA（1，1）模型分別對(duì)漢口站和沙市站的剩余序列建模，由矩估計(jì)法［26］估計(jì)模型未知參數(shù)，并求出相關(guān)系數(shù)，具體結(jié)果如表8所示?？梢钥闯觯夯▓@口站年徑流剩余序列由AIC和BIC準(zhǔn)則判定的階數(shù)為1階，但與其相依成分之間的相關(guān)系數(shù)為0.0240，比rα=0.2172小，因此該序列無(wú)相依變異，階數(shù)也最終確定為0；漢口站序列的相關(guān)系數(shù)r=0.3418，在區(qū)間[rα，rβ]內(nèi)，故其相依變異程度為弱變異；沙市站序列的相關(guān)系數(shù)r=0.4022，在區(qū)間[rβ，0.6] 內(nèi)，因此其相依變異程度為中變異。

表8 實(shí)測(cè)序列相依變異強(qiáng)弱分級(jí)

圖6為各站點(diǎn)實(shí)測(cè)徑流序列與其相依成分的示意圖。圖6（b）與圖6（c）為分別使用ARMA（2，1）模型和ARMA（1，2）模型對(duì)漢口站月徑流序列進(jìn)行擬合得到的原始水文序列及其相依成分示意圖，其中ARMA（2，1）模型的相關(guān)系數(shù)為0.3418、ARMA（1，2）模型的相關(guān)系數(shù)為0.3343，從相關(guān)系數(shù)的大小初步判斷ARMA（2，1）模型優(yōu)于ARMA（1，2）。進(jìn)一步引入一個(gè)能夠驗(yàn)證水文模型模擬結(jié)果優(yōu)劣的評(píng)價(jià)指標(biāo)即納什效率系數(shù)（NSE）［38］，其表達(dá)式為：

圖6 各站點(diǎn)原始水文序列及其相依成分示意

式中：T為序列長(zhǎng)度；Q0t為實(shí)測(cè)序列的第t個(gè)值；Qmt為模擬序列的第t個(gè)值；表示實(shí)測(cè)序列的均值。E的取值范圍為(-∞，1]，E越接近1表示擬合效果越好，模型可信度越高；E接近0，表示模擬結(jié)果總體可信，但擬合過(guò)程存在一定的誤差；E遠(yuǎn)小于0，則模型是不可信的。經(jīng)計(jì)算，使用AR?MA（2，1）模型和ARMA（1，2）模型對(duì)漢口站月徑流序列進(jìn)行擬合所得的納什效率系數(shù)分別為E1=-0.6104 和E2=-41.5755，考慮到漢口站序列的相依變異程度為弱變異，而相依變異程度越高，擬合效果越好，即擬合誤差越小，所以漢口站的相依成分序列與原序列間的擬合效果相比于相依中、強(qiáng)、巨變異要差，此為擬合過(guò)程存在誤差的原因，從而使用ARMA（2，1）模型對(duì)漢口站月徑流序列進(jìn)行擬合所得的納什效率系數(shù)為E1=-0.6104 得到了很好的解釋。故根據(jù)擬合效果與漢口站月徑流序列的相依變異程度可以判定漢口站月徑流的剩余序列符合ARMA（2，1）模型，這為AIC和BIC函數(shù)最小值不唯一或階數(shù)p和q有多種組合方式時(shí)，提供了一種對(duì)ARMA（p，q）相依成分進(jìn)行定階的新思路。圖6直觀地顯示，實(shí)測(cè)徑流序列與其相依成分的擬合程度逐漸加強(qiáng)，說(shuō)明實(shí)測(cè)序列的相依變異等級(jí)越來(lái)越高，故該分級(jí)方法的適用性在實(shí)測(cè)水文序列的應(yīng)用中得到了驗(yàn)證。

受自然和人為等因素的影響，水文要素常常呈現(xiàn)出一定的相依特性，故從氣候和人為兩方面對(duì)各站徑流序列的相依變異做出物理成因分析［21］?；▓@口站位于黃河下游，它的年徑流序列無(wú)相依變異主要是其總體受氣候變化與人類(lèi)活動(dòng)的影響較?。?3］；相較于氣候變化，人類(lèi)活動(dòng)對(duì)長(zhǎng)江中游徑流變化的影響更加顯著［39］，沙市水文站地處湖北省荊州市，作為長(zhǎng)江中游荊江河段的關(guān)鍵控制站，它上游葛洲壩和三峽水利樞紐的運(yùn)行和建設(shè)是導(dǎo)致此站月徑流中變異的重要原因，其次還受到沿程護(hù)岸工程以及洲灘守護(hù)工程的影響［40-41］；漢口站地處湖北武漢，是對(duì)漢江入?yún)R長(zhǎng)江后的水情進(jìn)行控制和監(jiān)測(cè)的重要站點(diǎn)，雖然受到長(zhǎng)江上游水利工程的影響，但漢江的水利工程設(shè)施相對(duì)較少，匯入長(zhǎng)江后，水文相依變異程度減輕為弱變異［13］。

5 結(jié)論

為定量研究水文序列中相依性的強(qiáng)弱，本文提出了一種基于自回歸滑動(dòng)平均模型（ARMA）的水文序列相依變異分級(jí)方法，分析驗(yàn)證后得出結(jié)論如下：

（1）由分級(jí)原理可知，ARMA模型中的自回歸系數(shù)和滑動(dòng)平均系數(shù)的大小影響著相依成分與水文時(shí)間序列間相關(guān)系數(shù)r的取值，而滑動(dòng)平均系數(shù)和自回歸系數(shù)的參數(shù)解與樣本自相關(guān)系數(shù)存在函數(shù)關(guān)系，基于此建立了相關(guān)系數(shù)r與能夠表征水文序列相依程度的樣本自相關(guān)系數(shù)之間的聯(lián)系，即相關(guān)系數(shù)r可以表示樣本序列的相依關(guān)系并描述其相依程度。此外，本文所推導(dǎo)的相關(guān)系數(shù)r與ARMA（p，q）模型參數(shù)間的關(guān)系（即式（19））在以往的文獻(xiàn)中并未提及，該關(guān)系可為后續(xù)對(duì)自回歸滑動(dòng)平均模型與相關(guān)系數(shù)的研究提供參考。

（2）運(yùn)用統(tǒng)計(jì)試驗(yàn)證明了分級(jí)原理與方法對(duì)較低階數(shù)的ARMA模型的合理性。①公式計(jì)算和統(tǒng)計(jì)試驗(yàn)所得r值差別很小，證明了式（19）的合理性；②在ARMA（1，1）模型中，相關(guān)系數(shù)r和自相關(guān)系數(shù)|ρ1|均隨自回歸系數(shù)的絕對(duì)值 |φ1|和滑動(dòng)平均系數(shù)的絕對(duì)值 |θ1|的增大而增大；在ARMA（1，2）模型中，r與 |ρ1|呈正相關(guān)關(guān)系、與 |ρ2|呈負(fù)相關(guān)關(guān)系，說(shuō)明以相關(guān)系數(shù)作為指標(biāo)衡量序列的相依性強(qiáng)弱是合理的；③根據(jù)整體變量H的取值范圍，生成模擬水文序列，所得相關(guān)系數(shù)處于對(duì)應(yīng)變異程度的取值區(qū)間，進(jìn)一步說(shuō)明此方法合理可靠；④而本文所提的分級(jí)原理與方法對(duì)于更高階數(shù)的ARMA（p，q）模型是否適用尚需進(jìn)一步研究與探討。

（3）實(shí)測(cè)序列分析表明：黃河花園口站年徑流序列無(wú)相依變異、長(zhǎng)江漢口（武漢關(guān)）站月徑流序列為相依弱變異、長(zhǎng)江沙市（二郎磯）站月徑流序列為中等程度的相依變異。并結(jié)合物理成因從氣候變化和人類(lèi)活動(dòng)兩個(gè)方面很好地解釋了各站點(diǎn)實(shí)測(cè)水文序列相依性依次加強(qiáng)的原因，驗(yàn)證了所提分級(jí)方法的合理性及對(duì)實(shí)測(cè)水文資料的適用性。

（4）提供了一種具有ARMA（p，q）相依成分的水文序列定階的新思路：在AIC和BIC函數(shù)最小值不唯一或階數(shù)p和q有多種組合方式時(shí)，可以先根據(jù)原始水文序列的相依變異程度初步判定模型階數(shù)是否合適，再根據(jù)其與相依成分的擬合效率系數(shù)最終確定模型階數(shù)。