秦耀軍，周曉勇，楊亞賓，賈 超

（1.山東省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局第二水文地質(zhì)工程地質(zhì)大隊(duì)，山東 德州 253072；2.山東大學(xué) 土建與水利學(xué)院，濟(jì)南250061）

工程勘察

基于數(shù)據(jù)同化技術(shù)的地質(zhì)參數(shù)反演分析研究

秦耀軍1，周曉勇2，楊亞賓1，賈 超2

（1.山東省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局第二水文地質(zhì)工程地質(zhì)大隊(duì)，山東 德州 253072；2.山東大學(xué) 土建與水利學(xué)院，濟(jì)南250061）

通過(guò)數(shù)據(jù)同化技術(shù)用卡爾曼濾波算法對(duì)觀測(cè)井所得水位監(jiān)測(cè)資料進(jìn)行同化計(jì)算，在得到較為精確水位觀測(cè)值的基礎(chǔ)上對(duì)主要取水層的滲透系數(shù)進(jìn)行參數(shù)反演，并與前期室內(nèi)試驗(yàn)獲得的滲透系數(shù)對(duì)比，表明通過(guò)數(shù)據(jù)同化過(guò)濾之后進(jìn)行的參數(shù)反演更精確有效。

數(shù)據(jù)同化；卡爾曼濾波；滲透系數(shù)；參數(shù)反演

1 研究背景

近年來(lái)隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，用數(shù)值分析方法研究地下水運(yùn)動(dòng)規(guī)律已發(fā)展成為地下水研究領(lǐng)域中的重要手段之一，特別是在對(duì)復(fù)雜的地下水問(wèn)題進(jìn)行分析時(shí)，數(shù)值分析方法能節(jié)省大量的人力物力。但為使刻畫(huà)的模型能夠符合真實(shí)情況，使計(jì)算過(guò)程加快收斂、縮減計(jì)算量，計(jì)算結(jié)果更加可靠，則需要提供精確的模型參數(shù)［1-3］。

含水層的地質(zhì)構(gòu)造雖然是確定的，但含水層的性質(zhì)在空間和時(shí)間上往往表現(xiàn)出極大的差異性，地下水流動(dòng)發(fā)生的環(huán)境一般又比較復(fù)雜，通過(guò)稀少的鉆孔資料很難對(duì)含水層特征進(jìn)行準(zhǔn)確的描述，所以根據(jù)掌握的觀測(cè)資料利用反演分析的方法反演出含水層的特征對(duì)研究地下水流動(dòng)十分必要。

由于觀測(cè)技術(shù)發(fā)展，易獲得很多原位觀測(cè)資料，但這些觀測(cè)資料都帶有一定誤差，而且多不連續(xù)，所以用這些直接觀測(cè)的數(shù)據(jù)對(duì)含水層性質(zhì)進(jìn)行描述也不精確，如何有效利用這些誤差較大的觀測(cè)數(shù)據(jù)，從中得到更有價(jià)值的信息，來(lái)降低數(shù)值模擬的不確定性，使之可以用于地下水的計(jì)算和預(yù)測(cè)，是一個(gè)十分重要的研究課題。目前常用數(shù)據(jù)同化技術(shù)來(lái)調(diào)和多來(lái)源的觀測(cè)信息，并對(duì)觀測(cè)數(shù)據(jù)進(jìn)行修正，最后從中得到較為準(zhǔn)確的信息，滿足正確刻畫(huà)地下水含水系統(tǒng)的要求。在數(shù)據(jù)同化方法中卡爾曼濾波法是應(yīng)用較廣，研究比較成熟的一種方法，最早由美國(guó)學(xué)者在20世紀(jì)70年代提出并用于軍事領(lǐng)域，經(jīng)過(guò)40多年的發(fā)展，已廣泛應(yīng)用在大氣、海洋、石油等多個(gè)領(lǐng)域［4-5］。

2006年，陳燕、張東曉等［6］首次運(yùn)用集合卡爾曼濾波（EnKF）同化地下水水頭及滲透系數(shù)的觀測(cè)資料，對(duì)模型的滲透系數(shù)場(chǎng)進(jìn)行了估計(jì)，這是集合卡爾曼濾波第一次在地下水領(lǐng)域中得到應(yīng)用。2008年，Liu和Chen［7］根據(jù)實(shí)驗(yàn)場(chǎng)的地下水流場(chǎng)，運(yùn)用該技術(shù)同化水文地質(zhì)觀測(cè)資料，對(duì)地下水模型參數(shù)做出了良好的估計(jì)。2009年，Huang［8］運(yùn)用EnKF改善了地下水污染物運(yùn)移的預(yù)測(cè)結(jié)果，說(shuō)明了EnKF相比傳統(tǒng)方法的優(yōu)勢(shì)。2010年，南通超、吳吉春等［9-13］通過(guò)局域化函數(shù)對(duì)小集合估計(jì)的協(xié)方差進(jìn)行修正降低了取樣噪聲對(duì)協(xié)方差估計(jì)的干擾，提高了濾波精度，并將其運(yùn)用到理想的二維承壓水模型中，得到很好的效果。2013年，史良勝等［13-15］將集合卡爾曼濾波方法拓展至補(bǔ)給條件下潛水流動(dòng)的數(shù)據(jù)同化問(wèn)題中，并通過(guò)同化水位、水力傳導(dǎo)度和降雨補(bǔ)給等測(cè)量數(shù)據(jù)來(lái)更新模型狀態(tài)、反演模型參數(shù)，得出在有長(zhǎng)期水位動(dòng)態(tài)測(cè)量數(shù)據(jù)時(shí)，可通過(guò)水位觀測(cè)值有效地反演出水力傳導(dǎo)度和降雨入滲補(bǔ)給系數(shù)的結(jié)論。2015年，林琳、史良勝等［16］指出用確定性集合卡爾曼濾波方法能夠有效減小大尺度問(wèn)題的抽樣誤差，通過(guò)算例表明該法能夠緩解在小樣本條件下的系統(tǒng)方差快速衰減，并探討了該方法在非均質(zhì)介質(zhì)和大尺度問(wèn)題中的應(yīng)用效果。通過(guò)近些年的發(fā)展集合卡爾曼濾波在地下水參數(shù)反演方面能夠獲取較好的結(jié)果，因此在地下水流場(chǎng)研究中被廣泛應(yīng)用。

2 數(shù)據(jù)同化技術(shù)的基本原理

卡爾曼濾波理論是1960被美國(guó)學(xué)者Kalman［6］和Bucy提出，對(duì)具有高斯誤差統(tǒng)計(jì)特性數(shù)據(jù)進(jìn)行最優(yōu)順序同化的一種方法。從包含誤差的數(shù)據(jù)資料里對(duì)所求變量，利用最小二乘估計(jì)、最小方差估計(jì)等方法進(jìn)行最佳估計(jì)，其實(shí)質(zhì)是由量測(cè)值重構(gòu)系統(tǒng)的狀態(tài)向量，以“預(yù)測(cè)—實(shí)測(cè)—修正”順序進(jìn)行遞推，得到關(guān)于狀態(tài)量的最優(yōu)估計(jì)。首先根據(jù)參考場(chǎng)和觀測(cè)值的特征誤差分布，對(duì)隨機(jī)場(chǎng)和得到的觀測(cè)值進(jìn)行一定擾動(dòng)，并利用帶有擾動(dòng)的隨機(jī)場(chǎng)和觀測(cè)值進(jìn)行同化分析計(jì)算，得到同化值，然后用這一組估計(jì)值的差異作為統(tǒng)計(jì)樣本進(jìn)行誤差協(xié)方差的估計(jì)。對(duì)這組估計(jì)值經(jīng)過(guò)預(yù)報(bào)后，又得到一組預(yù)報(bào)值，再將這組預(yù)報(bào)值重新作為下步同化的隨機(jī)場(chǎng)繼續(xù)進(jìn)行計(jì)算。集合卡爾曼濾波與卡爾曼濾波非常相似，其最大不同在于前者應(yīng)用了集合的思想，其誤差協(xié)方差矩陣是通過(guò)統(tǒng)計(jì)方式進(jìn)行計(jì)算，不需要時(shí)間上的遞推，而卡爾曼濾波算法中計(jì)算誤差協(xié)方差矩陣需要顯式計(jì)算并在時(shí)間步上進(jìn)行遞推。

地下水流動(dòng)是一個(gè)動(dòng)態(tài)變化的過(guò)程，在不同的時(shí)間有不同的觀測(cè)資料（如水頭值）加入。當(dāng)有帶有誤差的觀測(cè)資料加入時(shí)，卡爾曼濾波就可以將觀測(cè)資料的有用信息吸收到同化系統(tǒng)中去，連續(xù)獲得時(shí)間順序上的最優(yōu)估計(jì)值，并且不需要重新考慮前面己經(jīng)處理過(guò)的數(shù)據(jù)，所以地下水滲流用卡爾曼濾波進(jìn)行同化處理，通過(guò)不斷優(yōu)化觀測(cè)資料，對(duì)同化系統(tǒng)進(jìn)行實(shí)時(shí)修改，從而獲得含有更小誤差的估計(jì)值。

卡爾曼濾波算法主要包括預(yù)測(cè)和更新兩步，其主要計(jì)算步驟由式（1）～式（4）給出。

式（1）為EnKF的預(yù)報(bào)計(jì)算過(guò)程，X為狀態(tài)向量，F(xiàn)為預(yù)報(bào)算子，e為誤差，該式表示第i-1到第i時(shí)間步長(zhǎng)上的預(yù)報(bào)過(guò)程。將前一步同化值作為下一步的初值處理得到預(yù)報(bào)值。

式（2）為觀測(cè)向量的遞推過(guò)程，H為觀測(cè)算子，d為觀測(cè)向量，加入隨機(jī)誤差e后則與理論值相吻合。

式（3）為卡爾曼增益系數(shù)矩陣k的遞推過(guò)程，P為狀態(tài)誤差協(xié)方差向量，R為觀測(cè)誤差協(xié)方差向量。

3 工程算例分析

3.1 工程概況

研究區(qū)地層屬于魯西北沖洪積平原、沖海積平原區(qū)，根據(jù)該區(qū)域含水層介質(zhì)的巖性特征，本區(qū)域的水存在于第四系和松散沉積巖類(lèi)的裂隙孔隙中。受影響的因素較多，地質(zhì)構(gòu)造、地形地貌、地層巖性及氣象水文等都會(huì)對(duì)地下水的運(yùn)動(dòng)造成一定影響。根據(jù)現(xiàn)有勘測(cè)資料，可劃分為5個(gè)水文地質(zhì)亞區(qū)，地質(zhì)環(huán)境具有從山前到平原再到濱海的變化特點(diǎn)，含水巖組由全淡水的單層水質(zhì)結(jié)構(gòu)區(qū)逐步過(guò)渡到表層淡水～咸水～淡水的3層水質(zhì)結(jié)構(gòu)區(qū)乃至到濱海的全咸水區(qū)特點(diǎn)。

本文根據(jù)研究?jī)?nèi)容以德州市市區(qū)為計(jì)算范圍，依據(jù)前期搜集的資料分析可知該區(qū)取水層主要集中在300～500m的第3含水層到第5含水層，為此本文選取該取水段為研究對(duì)象對(duì)此段地層的參數(shù)進(jìn)行反演分析研究。

3.2 觀測(cè)資料的數(shù)據(jù)同化分析

本節(jié)基于卡爾曼濾波在地下水?dāng)?shù)據(jù)同化計(jì)算中的應(yīng)用，借助卡爾曼濾波計(jì)算程序?qū)τ^測(cè)孔水頭觀測(cè)值進(jìn)行同化計(jì)算，以獲得較為精確的水頭值，然后通過(guò)地下水流場(chǎng)分析軟件MODFLOW進(jìn)行參數(shù)反演。

根據(jù)最近德城區(qū)地下水開(kāi)采資料顯示，德城區(qū)深層地下水的主要開(kāi)采層位為300～500m含水層，開(kāi)采井的數(shù)量為170眼，占開(kāi)采井總數(shù)的65%，實(shí)際開(kāi)采總量2231.2萬(wàn)m3/a，允許開(kāi)采資源量1750萬(wàn)m3/a，開(kāi)采程度123%，屬于強(qiáng)超采，根據(jù)統(tǒng)計(jì)資料把調(diào)查所得開(kāi)采井放到模型中，根據(jù)獲得的觀測(cè)資料對(duì)每一口取水井水頭觀測(cè)值用卡爾曼濾波進(jìn)行過(guò)濾，部分水頭過(guò)濾結(jié)果如圖1～圖3。

圖1 ZK15觀測(cè)井水頭同化過(guò)程

圖2 ZK20觀測(cè)井水頭同化過(guò)程

圖3 誤差分析過(guò)程

3.3 研究區(qū)滲透系數(shù)反演分析

根據(jù)研究?jī)?nèi)容以德州市市區(qū)為計(jì)算范圍，依據(jù)前期搜集的資料分析可知該區(qū)取水層主要集中在300～500m的第3含水到第5含水層，為此選取該取水段為研究對(duì)象建立模型并進(jìn)行網(wǎng)格剖分如圖4。

圖4 研究區(qū)地下水分析數(shù)值模型及網(wǎng)格

共獲得83454個(gè)單元250362個(gè)節(jié)點(diǎn)，然后在建立模型的基礎(chǔ)上對(duì)研究區(qū)的地質(zhì)參數(shù)本文主要對(duì)滲透系數(shù)進(jìn)行反演。把各個(gè)觀測(cè)井過(guò)濾后的水位觀測(cè)值代入到反演分析軟件中，然后對(duì)研究的主要取水層的參滲透系數(shù)進(jìn)行反演，第3含水層到第5含水層滲透系數(shù)反演結(jié)果如圖5～圖10。

圖5 研究區(qū)第3含水層滲透性系數(shù)分布

圖6 研究區(qū)第3含水層滲透系數(shù)大小頻率

圖7 研究區(qū)第4含水層滲透系數(shù)分布

圖8 研究區(qū)第4含水層滲透系數(shù)大小頻率

圖9 研究區(qū)第5含水層滲透系數(shù)分布

圖10 研究區(qū)第5含水層滲透系數(shù)大小頻率

3.4 計(jì)算結(jié)果分析

經(jīng)過(guò)反演可得到第3層含水層水平方向滲透系數(shù)均值1.44m/d，最大值1.96m/d，最小值0.210m/d；第4含水層的滲透系數(shù)0.976m/d，最大值1.49m/d，最小值0.103m/d；第5含水層的平均滲透系數(shù)0.773m/d，最大值1.48m/d，最小值0.100m/d。從滲透系數(shù)反演計(jì)算結(jié)果可知，研究范圍內(nèi)不同含水層及同一含水層的不同地段的滲透系數(shù)并不一致，有的甚至相差比較大。從計(jì)算結(jié)果云圖中可以發(fā)現(xiàn)監(jiān)測(cè)井較為密集的地區(qū)反演的滲透系數(shù)值變幅較小。根據(jù)前期工程勘測(cè)資料及室內(nèi)研究可知各地層的滲透系數(shù)如表1。

表1 部分試驗(yàn)土樣的滲透系數(shù)

反演結(jié)果與室內(nèi)試驗(yàn)所得滲透系數(shù)對(duì)比如圖11。

圖11 各地層滲透系數(shù)反演與試驗(yàn)結(jié)果對(duì)比

從圖11中可以看出，反演所得值與試驗(yàn)值變化趨勢(shì)相對(duì)一致，由于從現(xiàn)場(chǎng)鉆孔獲取地下300～500m土樣會(huì)引起土樣擾動(dòng)，因此室內(nèi)試驗(yàn)所得滲透系數(shù)偏大，排除此方面干擾，可得反演所得值較為精確可靠。

4 結(jié)語(yǔ)

（1）本文所研究含水層主要為地下300～500m層位，地質(zhì)參數(shù)難以確定，因此首先要對(duì)幾個(gè)主要含水層的滲透系數(shù)進(jìn)行反演。在參數(shù)反演過(guò)程中引進(jìn)卡爾曼濾波算法，并運(yùn)用該算法對(duì)觀測(cè)井所獲得的水頭值進(jìn)行同化過(guò)濾，并把同化后的水頭值代入地下水分析軟件，對(duì)所研究含水層的滲透系數(shù)進(jìn)行了參數(shù)反演。

（2）通過(guò)反演計(jì)算得到第3含水層到第5含水層的滲透系數(shù)，計(jì)算結(jié)果表明各個(gè)含水的滲透系數(shù)不相同，并且同一含水層不同地方也相差較大，其中，第3含水層到第5含水層水平方向滲透系數(shù)均值分別為 1.44，0.97，0.77m/d，最大值分別為1.96，1.49，1.48m/d，最小值分別為0.2，0.1，0.1m/d。

（3）最后通過(guò)反演所得300～500m各地層滲透系數(shù)值與室內(nèi)試驗(yàn)研究所得滲透進(jìn)行對(duì)比分析，發(fā)現(xiàn)數(shù)據(jù)同化后再進(jìn)行反演所得值比較精確可靠。

［1］邵景力，趙宗壯，崔亞莉，等.華北平原地下水流模擬及地下水資源評(píng)價(jià)［J］.資源科學(xué)，2009，31（3）：361-367.

［2］謝先紅.基于數(shù)據(jù)同化方法的水文地質(zhì)參數(shù)反演與變量估計(jì)［A］.北京力學(xué)會(huì)第十六屆學(xué)術(shù)年會(huì)論文集［C］.2010.

［3］陳彥，吳吉春.含水層滲透系數(shù)空間變異性對(duì)地下水?dāng)?shù)值模擬的影響［J］.水科學(xué)進(jìn)展，2005（4）：482-487.

［4］施小清，吳吉春，袁永生.滲透系數(shù)空間變異性研究［J］.水科學(xué)進(jìn)展，2005（2）：210-215.

［5］Kalman R E.A new approach to linear filtering and prediction problems［J］.Transactions of the ASME-Journal of Basic Engineering，1960，82（series D）：35-45.

［6］Yan Chen， Dongxiao Zhang.Data assimilation for transient flow in geologic formations via ensemble Kalman filter［J］.Advances in Water Resources， 2006，29（8）：1107-1122.

［7］Liu C S， Chen Y， Zhang D X.Ivestigate of flow and transport processes at the MADE site using ensemble Kalman filter［J］.Advances in water resources，2008，31（7）：975-986.

［8］Huang C L， Hu B， Li X， et al.Using data assimilation method to calibrate a heterogeneous conductivity field and improve solute transportprediction with an unknown contamination source［J］.Stochastic Environmet Research and risk Assessment，2009，23（8）：1155-1167.

［9］Evensen G.The Ensemble Kalman Filter： theoretical f ormulation and practical implementation［J］.Ocean Dynamics，2003， 53（4）：1616-7228.

［10］Burgers G， Jan Van Leewen P， Evensen G.Analysis scheme in the ensemble kalman filter ［J］.Monthly weather review，1998，126（6）：1719-1724.

［11］南統(tǒng)超，吳吉春.集合卡爾曼濾波估計(jì)水文地質(zhì)參數(shù)的局域化修正［J］.水科學(xué)進(jìn)展，2010（5）：613-621.

［12］崔凱鵬，吳吉春.觀測(cè)數(shù)據(jù)時(shí)空密度對(duì)集合卡爾曼濾波計(jì)算精度的影響［J］.水利學(xué)報(bào)，2013（8）：915-923.

［13］崔凱鵬.集合卡爾曼濾波在地下水流及溶質(zhì)運(yùn)移數(shù)據(jù)同化中的應(yīng)用探討［D］.南京：南京大學(xué)，2013.

［14］史良勝，楊金忠，李少龍，等.基于KL-Galerkin解法的地下水流動(dòng)隨機(jī)分析［J］.四川大學(xué)學(xué)報(bào)（工程科學(xué)版），2005（5）：31-35.

［15］宋雪航，史良勝，楊金忠.基于集合卡爾曼濾波的潛水動(dòng)態(tài)預(yù)測(cè)方法［J］.武漢大學(xué)學(xué)報(bào)（工學(xué)版），2014（3）：324-331.

［16］林琳，史良勝，宋雪航.地下水參數(shù)反演的確定性集合卡爾曼濾波方法［J］.武漢大學(xué)學(xué)報(bào)（工學(xué)版），2016（2）：161-167，172.

Back analysis of geological parameters based on data assimilation

QIN Yao-jun1， ZHOU Xiao-yong2， YANG Ya-bin1， JIA Chao2
（1.NO.2 Hydro-engineering Geology Brigade of Shandong Provincial Bureau of Geology， Dezhou 253072， China；2.School of Civil Engineering， Shandong University， Jinan 250061， China）

Monitoring data of water level of observation wells are calculated by data assimilation technology and Calman filtering algorithm.With this method，the permeability coefficient of the main aquifer is inversed base on the accurate water level values，and compared with the permeability coefficient obtained by test.The results show that the inversion parameter after data assimilation and filtering is more accurate and effective.

data assimilation； Kalman filtering； permeability coefficient； parameter inversion

P64 文獻(xiàn)標(biāo)識(shí)碼：B 文章編號(hào)：1672-9900（2017）06-0078-05

2017-08-22

秦耀軍（1964-），男（漢族），山東濟(jì)寧人，高級(jí)工程師，主要從事水文地質(zhì)、工程地質(zhì)、環(huán)境地質(zhì)、地災(zāi)防治研究工作，（Tel）18561188280。

王艷肖）