鄒 曄，胡 瑩，朱世芳

（山東省魯南地質(zhì)工程勘察院，山東兗州272100）

地下水位動態(tài)序列分析與預(yù)報研究

鄒 曄*，胡 瑩，朱世芳

（山東省魯南地質(zhì)工程勘察院，山東兗州272100）

根據(jù)地下水水位動態(tài)資料的特點(diǎn)，將地下水位非平穩(wěn)時間序列分解為趨勢項(xiàng)函數(shù)、周期項(xiàng)函數(shù)和隨機(jī)項(xiàng)3部分，并疊加后建立組合數(shù)學(xué)模型，對地下水水位動態(tài)進(jìn)行預(yù)測。最后通過實(shí)例檢驗(yàn)，此種建模和預(yù)報理論取得了較好的效果。

地下水水位動態(tài)；非平穩(wěn)時間序列；預(yù)測

1 概述

地下水動態(tài)取決于地下水的補(bǔ)給、徑流與排泄條件，是受地形地貌、地層巖性、氣象和水文等自然因素和人工開采、灌溉、排水等人為因素綜合作用的結(jié)果。地下水動態(tài)反映了地下水要素隨時間變化的狀況，為了合理利用地下水或有效防范其危害，必須掌握地下水動態(tài)。

目前國內(nèi)外已經(jīng)有很多學(xué)者研究了采用隨機(jī)性數(shù)學(xué)模型進(jìn)行地下水水位動態(tài)的預(yù)測，其中受到大家公認(rèn)的方法主要包括回歸分析法、頻譜分析法和時間序列分析法等。在各種隨機(jī)模型中，無論是頻譜分析法還是時間序列分析法，都要求地下水動態(tài)序列滿足平穩(wěn)性條件，而實(shí)際中的地下水動態(tài)序列，由于人為干擾因素越來越大以及觀測時間的有限性，絕大部分序列都不滿足平穩(wěn)性條件，而是存在一個趨勢項(xiàng)。此外，由于各種隨機(jī)因素的作用，單一的周期加趨勢項(xiàng)模型必然會產(chǎn)生大量殘差，即隨機(jī)項(xiàng)，從而使預(yù)報的精度大大降低。為此本文嘗試運(yùn)用參數(shù)模型方法研究非平穩(wěn)性地下水動態(tài)序列。

2 非平穩(wěn)時間序列分析

2.1 非平穩(wěn)時間序列模型

所謂非平穩(wěn)時間系列，即表示統(tǒng)計性質(zhì)隨時間而異的存在非常廣泛的那一類物理數(shù)據(jù)，在實(shí)際問題中經(jīng)常遇到的時間序列一般都是非平穩(wěn)的。由于非平穩(wěn)時間序列的一般性和復(fù)雜性特點(diǎn)，對它至今尚無統(tǒng)一處理的一般方法。

設(shè)D（t）是時間t的確定性函數(shù)，R（t）是一個各態(tài)歷經(jīng)的平穩(wěn)隨機(jī)過程，通?？梢圆捎眉臃Ｐ蛠硌芯糠瞧椒€(wěn)過程，作為實(shí)際物理過程的近似，

由H（t）的測量數(shù)據(jù)，用一種統(tǒng)計的方法估計函數(shù)D（t）所含的一些參數(shù)，識別、提取、預(yù)報趨勢函數(shù)項(xiàng)D（t）。人們把這樣一類方法，稱為參數(shù)模型方法。

在時間序列分析中，識別、提取趨勢項(xiàng)D（t），是很重要的一項(xiàng)工作。為了提取D（t），一般假定D（t）由2部分組成，即：

式中：T（t）——實(shí)際測量數(shù)據(jù)隨時間t變化的主值函數(shù)項(xiàng)；

P（t）——由實(shí)測數(shù)據(jù)中可分離出來的周期函數(shù)項(xiàng)。

因此，一個非平穩(wěn)時間序列可以認(rèn)為是由趨勢成分、近似周期成分和平穩(wěn)隨機(jī)成分組成，而平穩(wěn)時間序列則要求統(tǒng)計參數(shù)的期望值與方差不隨時間改變。

2.2 非平穩(wěn)時間序列建模

2.2.1 建模思路

非平穩(wěn)水位動態(tài)序列H（t）是由趨勢成分T（t）、近似周期成分P（t）和平穩(wěn)隨機(jī)成分R（t）組成，其可表示為3個組成部分之和，表達(dá)式為[1]：

H（t）=T（t）+P（t）+R（t）

在非平穩(wěn)水位動態(tài)模型中，趨勢項(xiàng)T（t）反映變量的多年變化趨勢；周期項(xiàng)T（t）反映變量的周期性變化。趨勢項(xiàng)、周期項(xiàng)這2項(xiàng)反映了時間序列變量變化中的確定性成分，把這2項(xiàng)分離出去，余下的就是隨機(jī)項(xiàng)了。隨機(jī)項(xiàng)可用平穩(wěn)時間序列來分析，其可分為2項(xiàng)：平穩(wěn)時間序列項(xiàng)S（t）和純隨機(jī)項(xiàng)N（t），即：

R（t）=R（t）+N（t）

其中純隨機(jī)項(xiàng)N（t）作為白噪聲來處理。

將時間序列變量分解成3個組成部分之后，就可按各組成項(xiàng)的變化規(guī)律對未來時刻進(jìn)行外推，再將各項(xiàng)合成作出預(yù)報。

2.2.2 地下水水位動態(tài)時間序列模型的建立

（1）趨勢項(xiàng)分析：趨勢項(xiàng)是指在時間序列中，序列穩(wěn)定而規(guī)則的變動，即隨時間的推移對平均值來說增大或減小的趨勢，其可能是線性的，也可能是非線性的。趨勢項(xiàng)數(shù)學(xué)模型的結(jié)構(gòu)基本上由人們的經(jīng)驗(yàn)判定，當(dāng)無法判定趨勢項(xiàng)應(yīng)采用何種形式的數(shù)學(xué)結(jié)構(gòu)時，通常用確定性模型擬和趨勢成分T（t），用逐步回歸方法對模型系數(shù)加以求解。

對于趨勢分量T（t）可用多項(xiàng)式逼近，即：

由此，可采用多元回歸方法確定待定系數(shù)c0，c1，c2，…，c10和階數(shù)。其具體求解方法通過編程序統(tǒng)一處理來實(shí)現(xiàn)。如果經(jīng)逐步回歸計算，回歸系數(shù)全為零，可以認(rèn)為（3）式無趨勢項(xiàng)T（t）。

（2）周期項(xiàng)分析：

①頻譜分析方法：周期分量是序列隨著時間的推移而呈現(xiàn)出的周期性成分。時間序列分離趨勢之后，將剩余的序列P（t）=H（t）-T（t）進(jìn)行周期分析。識別和提取周期項(xiàng)P（t）的方法有方差分析、頻譜分析和周期圖分析等，這里主要采用頻譜分析法。

頻譜分析是利用傅立葉級數(shù)把某個資料的時間序列表示成無數(shù)個不同周期的簡諧波和的形式來分析序列變化規(guī)律的一種方法。對序列P（t）可用L個波疊加的形式表示其周期項(xiàng)[2]：

其中的每個項(xiàng)為一個分波，分別稱Ai，ωi，φi為第i個分波的振幅、頻率和相位，a0為一常數(shù)。因?yàn)?/p>

sin( ) ωit+φi=sinωitcosφi+cosωitsinφi

并令：

ai=Aisinφibi=Aicosφi

則P（t）可以改記為：

由于觀測資料的有限性，無法進(jìn)行無窮分波，因此一般假定P（t）有K個分波（試驗(yàn)周期個數(shù)），即：

對于給定的時間序列可用以下的方法來確定a0，ai，bi，i=1，2，…，K。設(shè)有n個水位H（t），t=1，2，…，n，n為樣本的長度，除去其趨勢項(xiàng)T（t）后，余下的序列記為P（t），即P（t）=H（t）-T（t），t=1，2，…，n，又認(rèn)為K個分波各有年的周期，即第i個分波，周期為：

因而第i個分波的頻率為：

于是：

可以采用最小二乘法來確定系數(shù)a0，ai，bi，通過一系列推導(dǎo)，可以求得：

如何從上述求得的各分波系數(shù)ai，bi中選取主要周期。常用下面2種方法推求。

②主要周期判斷：本次分析用振幅的大小判斷主要周期。

因：

故：

在顯著水平為0.05時，若：

則認(rèn)為相應(yīng)的第i個分波是主要周期，上式中的n為樣本長度，K為試驗(yàn)周期個數(shù)，σ2為序列P（t）的方差，可用樣本方差S2來估計。通過推導(dǎo)能夠證明S2可用下式來計算：

其中：

若A1，…，Ak中有M個分波滿足（9）式，則我們在（6）式中只保留相應(yīng)的M個項(xiàng)即可，即最后求得P（t）的周期項(xiàng)為：

在實(shí)際應(yīng)用中為節(jié)省工作量，通常在L個波中選取波動比較顯著的幾個諧波相加得到周期項(xiàng)，一般只需選取前6個顯著諧波就能滿足精度要求了。

（3）隨機(jī)項(xiàng)分析：隨機(jī)分量是指各時刻的序列值與前一個時刻或幾個時刻值存在著某種相關(guān)關(guān)系的成分。時間序列模型中，消除趨勢項(xiàng)和近似周期項(xiàng)后的剩余序列，即為平穩(wěn)隨機(jī)系列項(xiàng)，此時，即可用平穩(wěn)隨機(jī)模型方法來分析求解。

①求解數(shù)學(xué)模型：在數(shù)理統(tǒng)計中，我們已知回歸模型為：

它表示觀測值yt對另一組觀測值（x1t，x2t，…，xkt）的相依性，上式可以視為由2部分組成，一部分取決于自變量（x1t，x2t，…，xkt），而另一部分則是隨機(jī)成分εt。觀測序列｛yt｝是假定為相互獨(dú)立或不相關(guān)的。

對時間序列xt=x（t）（假定屬平穩(wěn)序列，且均值為零）有[3]：

上式為p階自回歸模型，常記為AR（p），式中φ1，φ2，…，φp為待定常數(shù)，稱它們?yōu)樽曰貧w系數(shù)，εt表示在量測過程中存在的隨機(jī)干擾和未來預(yù)報中出現(xiàn)的誤差。我們假定t時的觀測誤差εt與t以前的觀測值xt-j，j=1，2，…，p獨(dú)立，因而有：

②自回歸系數(shù)的確定：自回歸系數(shù)φ1，φ2，…，φp的確定有2條途徑。

其一是用xt-i乘（14）式兩邊，得到：

上式兩邊取數(shù)學(xué)期望，并利用（15）式得：

此即一組p階線性方程：

式中：R（1），R（2），…，R（p）可以用估計值r（1），r（2），…，r（p）來代替，于是從理論上講由（18）式便可以解出φ1，φ2，…，φp。

根據(jù)x1，x2，…，xN為各態(tài)歷經(jīng)平穩(wěn)時間序列的假定，可以用下式估計相關(guān)函數(shù)r（τ）[4]：

若x1，x2，…，xp未標(biāo)準(zhǔn)化，需先標(biāo)準(zhǔn)化如下：

然后對xb1，xb2，…，xbN使用（19）式。這樣，我們得到自回歸系數(shù)φ1，φ2，…，φp的相應(yīng)估計值b1，b2，…，bp，它們滿足p階方程組：

式中r（p）[r（τ）]即為相關(guān)函數(shù)R（p）[R（τ）]的估計值，見式（21）。

另一途徑用最小二乘估計的方法來決定（14）式中的系數(shù)φ1，φ2，…，φp，即選φ1，φ2，…，φp，使：

為最小值[5]。上式關(guān)于φj，j=1，2，…，p分別求偏導(dǎo)數(shù)，得：

由各態(tài)歷經(jīng)性的假定有：

于是（22）式兩邊除以N-p，據(jù)（23）式，便得：

將φk改記為bk再展開后，它與（20）式完全一樣。

用各種代數(shù)方法求解方程組，得φ1，φ2，…，φp。

③自回歸模型階數(shù)p的確定：上面給出了自回歸模型（14）式，從上面的過程可以看出，如果取的階數(shù)p不同，則所得的模型也不同，因此便提出一個問題：究竟p應(yīng)是多少，才可使得求出的自回歸模型（14）最佳？這實(shí)際上是個模型識別問題。

在使用上可以這樣考慮：給出一個誤差范圍，例如0.01，即（14）中的 ||εt＜0.01，然后對不同的p進(jìn)行計算，直至滿足 ||εt＜0.01的那個p為止，最后所得結(jié)果便是所要的值。

3 實(shí)例驗(yàn)證

為了驗(yàn)證文中所提理論，筆者收集了某鉆孔1996～2012年的1224個水位值（每個月6個水位觀測值）作為建模數(shù)據(jù)，利用文章所提理論確定了模型，趨勢分析模型參數(shù)C0取3.510，C1取0.0042，其它模型參數(shù)結(jié)果見表1和表2。模型對1994～2012年水位進(jìn)行了擬合，取得了較好的效果，實(shí)際水位與擬合水位對比圖見圖1，同時利用建立的模型對該鉆孔2013～2027年的水位進(jìn)行了預(yù)報（見圖1）。

圖1 地下水位埋深動態(tài)擬合、驗(yàn)證、預(yù)測曲線

表1 周期分析參數(shù)表

表2 自回歸模型參數(shù)表

4 結(jié)論

非平穩(wěn)時間系列方法適宜樣本容量較大的地下水動態(tài)建模和預(yù)測，以容量大于100以上為佳，小樣本不足以全面反映地下水位動態(tài)的變化規(guī)律，不能有效地提取出真實(shí)的趨勢項(xiàng)和周期項(xiàng)，且受隨機(jī)干擾較大。同時該方法雖然在一定程度上考慮到了人為干擾的影響，但在偶然、大強(qiáng)度干擾下，預(yù)報結(jié)果會有一些誤差。

[1]張偉,徐建華,秦勇.非平穩(wěn)性地下水動態(tài)序列分析及預(yù)測[J].工程勘察，2000(1):17.

[2]燕良東.地下水動態(tài)組合模型研究[D].遼寧：遼寧工程技術(shù)大學(xué)，2011.

[3]丁晶,劉全授.隨機(jī)水文學(xué)[M].北京：中國水利水電出版社，1997.

[4]安鴻志,時間序列分析[M].上海：華東師范大學(xué)出版社，1992.

[5]申鼎煊，隨機(jī)過程[M].湖北：華中理工大學(xué)出版社，1990.

TV641.74

A

1004-5716(2015)06-0182-04

2014-06-25

2014-07-03

鄒曄（1982-），男（漢族），江蘇宜興人，工程師，現(xiàn)從事水、工、環(huán)技術(shù)工作。