孫永泉,劉梅俠,董學(xué)良,陳興國(guó)

(1.黑龍江省904水文地質(zhì)工程地質(zhì)勘察院,雙城 150100;2.大慶市環(huán)宇地質(zhì)工程有限公司,黑龍江 大慶 163000)

1 Kalman濾波技術(shù)的基本原理

1.1 Kalman濾波技術(shù)的概念

Kalman濾波技術(shù)[1]是指用觀測(cè)數(shù)據(jù),采用線形系統(tǒng)狀態(tài)方程和測(cè)量方程作出某一時(shí)刻的測(cè)量值的最優(yōu)估計(jì),并希望估計(jì)具有遞推解。

1.2 Kalman濾波技術(shù)的基本原理

卡爾曼濾波有兩個(gè)基本方程[2]:即狀態(tài)方程和測(cè)量方程。對(duì)地下水流系統(tǒng)而言。狀態(tài)方程是描述地下水流系統(tǒng)確定—隨機(jī)性問(wèn)題的方程式,若定義地下水流系統(tǒng)變量為水頭矢量H,則地下水流系統(tǒng)的狀態(tài)方程式為:

地下水流系統(tǒng)的測(cè)量方程為:

式中Yk為k時(shí)刻地下水系統(tǒng)狀態(tài)的測(cè)量矢量;Vk為k時(shí)刻地下水系統(tǒng)狀態(tài)的測(cè)量誤差矢量;為k時(shí)刻地下水系統(tǒng)狀態(tài)的測(cè)量矩陣,由監(jiān)測(cè)孔相對(duì)于模型分節(jié)點(diǎn)的位置確定。

式(1)和式(2)構(gòu)成地下水流系統(tǒng)確定—隨機(jī)性數(shù)值模型。該模型的求解可按卡爾曼濾波算法求得,其基本假設(shè)條件為[3]:

1)系統(tǒng)噪聲向量的數(shù)學(xué)期望為已知函數(shù),即:

2)測(cè)量誤差向量的期望值為零,即

3)測(cè)量誤差的協(xié)方差陣Rk為

4)系統(tǒng)噪聲的協(xié)方差陣Qk為

5)系統(tǒng)噪聲與測(cè)量誤差向量相互獨(dú)立,即:

在上式假設(shè)條件下,經(jīng)過(guò)一系列數(shù)學(xué)推導(dǎo),得到卡爾曼濾波最優(yōu)估計(jì)的遞推算法,即:

及地下水流系統(tǒng)模擬遞推算法為:

1.3 Kalman濾波技術(shù)應(yīng)用于地下水流系統(tǒng)確定—隨機(jī)性模型校正

運(yùn)用卡爾曼濾波最優(yōu)估計(jì)算法,對(duì)已建立的模型進(jìn)行校正,以確定系統(tǒng)參數(shù)A、B及Q。對(duì)于有限元法的不規(guī)則三角剖分而言[4]:

式中G與D分別為有限元方程中的導(dǎo)水矩陣及貯水矩陣;B1為輸入項(xiàng)有關(guān)的矩陣。通過(guò)求解A、B可確定地下水流系統(tǒng)的確定性參數(shù)K(或T),S(或μ)等。

式中ρ l為兩個(gè)節(jié)點(diǎn)間(i,j)單位距離上的系統(tǒng)噪聲相關(guān)系數(shù);dij為i,j兩點(diǎn)之間的距離,當(dāng)i=j時(shí),dij=0。

在進(jìn)行模型校正過(guò)程中,一般按正演方法求解[5],事先給定K,S,ρ1,σ2的初值,運(yùn)用卡爾曼濾波最優(yōu)估計(jì)算法計(jì)算,可得到新息矢量Nk、新息的理論方差矢量Z2、樣本方差矢量S2及樣本均值矢量ˉN,可根據(jù)以下統(tǒng)計(jì)分析方法,在給定參數(shù)估計(jì)的可靠度后,按概率統(tǒng)計(jì)方法判別模型校正是否滿意[6-7]。

由于新息矢量可表達(dá)為:

從理論上講,新息矢量的數(shù)學(xué)期望應(yīng)滿足:

新息矢量的理論方差為:

新息矢量的樣本均值為:

新息矢量的樣本方差為:

理論上講,若給定含水層系統(tǒng)確定性參數(shù)(K,S等),隨機(jī)性參數(shù)能使新息的均值和方差滿足下式:

則所確定的參數(shù)是最佳模型參數(shù),參數(shù)的可靠性為100%,此時(shí)模型校正結(jié)束。實(shí)際中由于樣本數(shù)NT有限及計(jì)算誤差及系統(tǒng)誤差的存在,實(shí)際上求得ˉN≠0,S2≠Z2。為了給 ˉN,S2賦予一個(gè)明確的精度概念,引入均值和方差的置信區(qū)間,只要 ˉN,S2在一定的置信區(qū)間內(nèi),則認(rèn)為模型校正結(jié)束,所確定的參數(shù)滿足一定的精度。

1.3.1 用S2確定總體均值的置信區(qū)間

新息矢量Nk為總體N(0,Z2)的一簇樣本,用樣本方差S2估計(jì)Z2,則[8-9]:

且t(NT-1)分布不依賴于總體均值,從而得:

于是,得均值為零(Nk=0)的置信度為1-α的置信區(qū)間為:

1.3.2 方差Z2的置信區(qū)間

由于

且χ2(NT-1)分布與Z2無(wú)關(guān),則有:

于是,得方差Z2置信度為1-α置信區(qū)間為:

如果全部均值和方差都通過(guò)檢驗(yàn),則可得出結(jié)論,卡爾曼濾波與地下水流數(shù)值耦合模型已得到很好地校正[10],其可信度為1-α。

2 實(shí) 例

2.1 地下水流系統(tǒng)數(shù)學(xué)模型

研究區(qū)面積約8 000 km2,主要取水層位為第三系承壓含水層。歷經(jīng)30 a的開采,已形成了4 000 km2的地下水位降落漏斗,并積累了大量的動(dòng)態(tài)資料。研究區(qū)東部為第三系承壓水缺失邊界,南部、西部和北部為已知水頭邊界。研究區(qū)地下水徑流方向?yàn)楸北蔽飨?主要補(bǔ)給方式為側(cè)向徑流補(bǔ)給,其次為垂向補(bǔ)給,主要排泄方式為人工開采,依據(jù)上述含水層系統(tǒng)特征,研究區(qū)滲流問(wèn)題的數(shù)學(xué)模型可描述為[11]:

式中Wh(x,y,t)為地下水流模型中的隨機(jī)誤差項(xiàng),也稱系統(tǒng)噪聲,它包括實(shí)際地下水系統(tǒng)概化成模型引起的誤差、邊界條件概化、初始條件概化引起的誤差等項(xiàng);H(x,y,t)為滲流域中任意時(shí)刻的地下水頭分布值(或稱地下水流系統(tǒng)的狀態(tài))(m);K M為含水層的導(dǎo)水系數(shù)(m2/d),K為含水層的滲透系數(shù)(m/d),M為含水層的厚度(m);S為含水層的貯水系數(shù),無(wú)量綱;H0(x,y)為滲流域中初始水頭分布或稱地下水流系統(tǒng)的初始狀態(tài)(m);Hb(x,y,t)為滲流域中已知邊界水頭分布(m);Qi為滲流域中i節(jié)點(diǎn)上源匯項(xiàng);Ω為滲流域;Γ1為滲流域中已知水頭(第一類)邊界;Γ2為滲流域中已知流量(第二類)邊界。

研究區(qū)離散為688個(gè)三角形單元,383個(gè)節(jié)點(diǎn),初步劃分為18個(gè)確定性參數(shù)區(qū)和13個(gè)隨機(jī)性參數(shù)區(qū)。

2.2 模型校正

運(yùn)用前述的卡爾曼濾波最優(yōu)估計(jì)遞推算法,選取已知監(jiān)測(cè)孔20個(gè)(選取原則:①保證每個(gè)參數(shù)區(qū)至少有1個(gè)監(jiān)測(cè)孔;②實(shí)測(cè)值中缺失較多者不選),計(jì)算時(shí)段選在1986年1月至1993年12月,在此時(shí)間段內(nèi)監(jiān)測(cè)孔數(shù)量較多,監(jiān)測(cè)序列連續(xù)性好。測(cè)量誤差的方差取為0.002 5 m2。校正后的滲透系數(shù)、貯水系數(shù)和系統(tǒng)噪聲標(biāo)準(zhǔn)差見(jiàn)表1,表2,系數(shù)ρl的校正值為0.85。

表1 校正后的模型確定性參數(shù)表T able 1 Certaintyof the model parameter table afer correcting

表2 校正后的模型的隨機(jī)性參數(shù)表T able 2 Randomness of the model parameter table afer correcting

圖1、圖2分別為研究區(qū)初始流場(chǎng)等值線圖和估計(jì)方差等值線圖及最優(yōu)估計(jì)流場(chǎng)和標(biāo)準(zhǔn)差等值線圖,從圖2a和圖2b(為最優(yōu)估計(jì)誤差標(biāo)準(zhǔn)差圖)可以看出,計(jì)算的地下水流場(chǎng)和實(shí)際的地下水流場(chǎng)擬合的比較好,因此,計(jì)算的地下水流場(chǎng)基本上代表了地下水位的空間分布特征,說(shuō)明模型正確地描述了含水層系統(tǒng)的基本特征。同時(shí)對(duì)20個(gè)監(jiān)測(cè)孔的歷時(shí)曲線進(jìn)行了對(duì)比,除7119、7101、7208 3個(gè)孔擬合較差外,其余17個(gè)孔擬合程度很好。每個(gè)監(jiān)測(cè)孔新息的均值見(jiàn)表3,表4和表5為模型校正可靠性的檢驗(yàn)結(jié)果。

由表3可見(jiàn):7119、7101、7208、8208、8409 5個(gè)孔的新息值較大。由表 4的 t檢驗(yàn)可見(jiàn):7119、7101、7208、8208及8409 5個(gè)孔未通過(guò)檢驗(yàn),可能是由于觀測(cè)精度及觀測(cè)序列中有缺失數(shù)據(jù)引起的,其余15個(gè)孔新息的均值與零無(wú)明顯差異。由表5可見(jiàn),除7114、7207、8208 3個(gè)孔未通過(guò)χ2檢驗(yàn)外,其余17個(gè)孔全部通過(guò)χ2檢驗(yàn),這幾個(gè)孔未通過(guò)χ2檢驗(yàn)的主要原因是系統(tǒng)噪聲標(biāo)準(zhǔn)差高及新息樣本標(biāo)準(zhǔn)差太大。要想全部通過(guò)t和χ2檢驗(yàn)是十分困難的,只要大多數(shù)孔通過(guò)檢驗(yàn)就可以了。通過(guò)模型校正,可認(rèn)為校正后的模型基本上反映了實(shí)際的水文地質(zhì)特征及地下水流系統(tǒng)特征,可靠性達(dá)90%以上(模型正檢驗(yàn)時(shí)選取置信度95%)。

表3 已知監(jiān)測(cè)孔新息的均值表Table 3 Known monitoring holes mean of the new interest rate table /m

表4 已知監(jiān)測(cè)孔的新息均值的t檢驗(yàn)T able 4 Known monitoring holes of the newrate means-t test

表5 已知監(jiān)測(cè)孔的新息標(biāo)準(zhǔn)差的檢驗(yàn)Table 5 Known monitoring holes test of the new interest rate standard deviation

3 結(jié) 論

本文運(yùn)用了卡爾曼濾波技術(shù)校正地下水流系統(tǒng),提出一種估計(jì)地下水流系統(tǒng)的隨機(jī)性參數(shù)的卡爾曼濾波耦合算法。通過(guò)對(duì)研究區(qū)地下水流系統(tǒng)的分析,可以得出以下結(jié)論:

1)本文的方法既考慮了地下水系統(tǒng)的確定性,又考慮了隨機(jī)性。

2)本文的方法可估計(jì)地下水流系統(tǒng)分布式參數(shù),包括確定性參數(shù)和隨機(jī)性參數(shù),并可給出參數(shù)估計(jì)的置信區(qū)間和精度。

3)估計(jì)的研究區(qū)地下水流系統(tǒng)分布型參數(shù)的可靠性達(dá)90%以上。

[1] 李官錫,賈洪利,周?？?等.大慶西部地區(qū)地下水補(bǔ)給條件及水質(zhì)的研究[R].大慶:大慶石油管理局供水公司,1991.

[2] 李俊亭.地下水流數(shù)值模擬[M].北京:地質(zhì)出版社,1989.

[3] J.Bear.地下水水力學(xué)[M].北京:地質(zhì)出版社,1985.

[4] 仵彥卿.地下水系統(tǒng)的基本概念與組成[J].長(zhǎng)安大學(xué)學(xué)報(bào)(地球科學(xué)版),1990,12(4):88-91.

[5] 戴冠中,佟明安.現(xiàn)代控制理論導(dǎo)論[M].北京:國(guó)防工業(yè)出版社,1989.

[6] 高鐘毓.工程系統(tǒng)中的隨機(jī)過(guò)程——隨機(jī)系統(tǒng)分析與最優(yōu)濾波[M].北京:清華大學(xué)出版社,1989.

[7] 浙江大學(xué)數(shù)學(xué)系高等數(shù)學(xué)教研組.概率論與數(shù)理統(tǒng)計(jì)[M].北京:高等教育出版社,1979.

[8] F.C.Van Geer,Van der Kloet P.Estimation of Parameters in Groundwater Flow Problems Using a Kalman Filter Algorithm [J].Symposium International,TNO Comm.Hydrol.Res.T he Hague,1983:449-462.

[9] F.C.Van Geer,C.B.M.Te Stroet,Zhou Yangxiao.Using Kalman Filtering to Improve and Quantify the U ncertainty of Numerical Groundwater Simulations 1:T he Role of System Noise and Its Calibration[J].Water Resour.Res.,1991,27 (8):1 987-1 994.

[10] F.C.Van Geer.Applicatiop of Kalman Filtering in the Analysis and Design of Groundwater Monitoring Networks[R].Report,Delft Univ.of T echnol.,Delft,Netherlands,1987.

[11] Zhou Yangx iao,C.B.M.T e Stroet,F.C.Van Geer.Using Kalman Filtering to Improve and Quantity the Uncertainty of Numerical Groundwater Simulation 2:Application to M onitoring Network Design[J].Water Resource.Res.,1991,27(8): 1 995-2 006.