鄒 曄，胡 瑩，朱世芳

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

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

鄒 曄*，胡 瑩，朱世芳

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

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

地下水水位動(dòng)態(tài)；非平穩(wěn)時(shí)間序列；預(yù)測(cè)

1 概述

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

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

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

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

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

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

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

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

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

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

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

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

2.2.1 建模思路

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

并令：

ai=Aisinφibi=Aicosφi

則P（t）可以改記為：

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

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

因而第i個(gè)分波的頻率為：

于是：

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

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

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

因：

故：

在顯著水平為0.05時(shí)，若：

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

其中：

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

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

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

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

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

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

上式為p階自回歸模型，常記為AR（p），式中φ1，φ2，…，φp為待定常數(shù)，稱它們?yōu)樽曰貧w系數(shù)，εt表示在量測(cè)過程中存在的隨機(jī)干擾和未來預(yù)報(bào)中出現(xiàn)的誤差。我們假定t時(shí)的觀測(cè)誤差εt與t以前的觀測(cè)值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）可以用估計(jì)值r（1），r（2），…，r（p）來代替，于是從理論上講由（18）式便可以解出φ1，φ2，…，φp。

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

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

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

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

另一途徑用最小二乘估計(jì)的方法來決定（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不同，則所得的模型也不同，因此便提出一個(gè)問題：究竟p應(yīng)是多少，才可使得求出的自回歸模型（14）最佳？這實(shí)際上是個(gè)模型識(shí)別問題。

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

[4]安鴻志,時(shí)間序列分析[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ù)工作。