基于數據同化技術的地質參數反演分析研究

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

（1.山東省地質礦產勘查開發局第二水文地質工程地質大隊，山東 德州 253072；2.山東大學 土建與水利學院，濟南250061）

工程勘察

基于數據同化技術的地質參數反演分析研究

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

（1.山東省地質礦產勘查開發局第二水文地質工程地質大隊，山東 德州 253072；2.山東大學 土建與水利學院，濟南250061）

通過數據同化技術用卡爾曼濾波算法對觀測井所得水位監測資料進行同化計算，在得到較為精確水位觀測值的基礎上對主要取水層的滲透系數進行參數反演，并與前期室內試驗獲得的滲透系數對比，表明通過數據同化過濾之后進行的參數反演更精確有效。

數據同化；卡爾曼濾波；滲透系數；參數反演

1 研究背景

近年來隨著計算機技術的快速發展，用數值分析方法研究地下水運動規律已發展成為地下水研究領域中的重要手段之一，特別是在對復雜的地下水問題進行分析時，數值分析方法能節省大量的人力物力。但為使刻畫的模型能夠符合真實情況，使計算過程加快收斂、縮減計算量，計算結果更加可靠，則需要提供精確的模型參數［1-3］。

含水層的地質構造雖然是確定的，但含水層的性質在空間和時間上往往表現出極大的差異性，地下水流動發生的環境一般又比較復雜，通過稀少的鉆孔資料很難對含水層特征進行準確的描述，所以根據掌握的觀測資料利用反演分析的方法反演出含水層的特征對研究地下水流動十分必要。

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

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

2 數據同化技術的基本原理

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

地下水流動是一個動態變化的過程，在不同的時間有不同的觀測資料（如水頭值）加入。當有帶有誤差的觀測資料加入時，卡爾曼濾波就可以將觀測資料的有用信息吸收到同化系統中去，連續獲得時間順序上的最優估計值，并且不需要重新考慮前面己經處理過的數據，所以地下水滲流用卡爾曼濾波進行同化處理，通過不斷優化觀測資料，對同化系統進行實時修改，從而獲得含有更小誤差的估計值。

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

式（1）為EnKF的預報計算過程，X為狀態向量，F為預報算子，e為誤差，該式表示第i-1到第i時間步長上的預報過程。將前一步同化值作為下一步的初值處理得到預報值。

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

式（3）為卡爾曼增益系數矩陣k的遞推過程，P為狀態誤差協方差向量，R為觀測誤差協方差向量。

3 工程算例分析

3.1 工程概況

研究區地層屬于魯西北沖洪積平原、沖海積平原區，根據該區域含水層介質的巖性特征，本區域的水存在于第四系和松散沉積巖類的裂隙孔隙中。受影響的因素較多，地質構造、地形地貌、地層巖性及氣象水文等都會對地下水的運動造成一定影響。根據現有勘測資料，可劃分為5個水文地質亞區，地質環境具有從山前到平原再到濱海的變化特點，含水巖組由全淡水的單層水質結構區逐步過渡到表層淡水～咸水～淡水的3層水質結構區乃至到濱海的全咸水區特點。

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

3.2 觀測資料的數據同化分析

本節基于卡爾曼濾波在地下水數據同化計算中的應用，借助卡爾曼濾波計算程序對觀測孔水頭觀測值進行同化計算，以獲得較為精確的水頭值，然后通過地下水流場分析軟件MODFLOW進行參數反演。

根據最近德城區地下水開采資料顯示，德城區深層地下水的主要開采層位為300～500m含水層，開采井的數量為170眼，占開采井總數的65%，實際開采總量2231.2萬m3/a，允許開采資源量1750萬m3/a，開采程度123%，屬于強超采，根據統計資料把調查所得開采井放到模型中，根據獲得的觀測資料對每一口取水井水頭觀測值用卡爾曼濾波進行過濾，部分水頭過濾結果如圖1～圖3。

圖1 ZK15觀測井水頭同化過程

圖2 ZK20觀測井水頭同化過程

圖3 誤差分析過程

3.3 研究區滲透系數反演分析

根據研究內容以德州市市區為計算范圍，依據前期搜集的資料分析可知該區取水層主要集中在300～500m的第3含水到第5含水層，為此選取該取水段為研究對象建立模型并進行網格剖分如圖4。

圖4 研究區地下水分析數值模型及網格

共獲得83454個單元250362個節點，然后在建立模型的基礎上對研究區的地質參數本文主要對滲透系數進行反演。把各個觀測井過濾后的水位觀測值代入到反演分析軟件中，然后對研究的主要取水層的參滲透系數進行反演，第3含水層到第5含水層滲透系數反演結果如圖5～圖10。

圖5 研究區第3含水層滲透性系數分布

圖6 研究區第3含水層滲透系數大小頻率

圖7 研究區第4含水層滲透系數分布

圖8 研究區第4含水層滲透系數大小頻率

圖9 研究區第5含水層滲透系數分布

圖10 研究區第5含水層滲透系數大小頻率

3.4 計算結果分析

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

表1 部分試驗土樣的滲透系數

反演結果與室內試驗所得滲透系數對比如圖11。

圖11 各地層滲透系數反演與試驗結果對比

從圖11中可以看出，反演所得值與試驗值變化趨勢相對一致，由于從現場鉆孔獲取地下300～500m土樣會引起土樣擾動，因此室內試驗所得滲透系數偏大，排除此方面干擾，可得反演所得值較為精確可靠。

4 結語

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

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

（3）最后通過反演所得300～500m各地層滲透系數值與室內試驗研究所得滲透進行對比分析，發現數據同化后再進行反演所得值比較精確可靠。

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

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

［3］陳彥，吳吉春.含水層滲透系數空間變異性對地下水數值模擬的影響［J］.水科學進展，2005（4）：482-487.

［4］施小清，吳吉春，袁永生.滲透系數空間變異性研究［J］.水科學進展，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］南統超，吳吉春.集合卡爾曼濾波估計水文地質參數的局域化修正［J］.水科學進展，2010（5）：613-621.

［12］崔凱鵬，吳吉春.觀測數據時空密度對集合卡爾曼濾波計算精度的影響［J］.水利學報，2013（8）：915-923.

［13］崔凱鵬.集合卡爾曼濾波在地下水流及溶質運移數據同化中的應用探討［D］.南京：南京大學，2013.

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

［15］宋雪航，史良勝，楊金忠.基于集合卡爾曼濾波的潛水動態預測方法［J］.武漢大學學報（工學版），2014（3）：324-331.

［16］林琳，史良勝，宋雪航.地下水參數反演的確定性集合卡爾曼濾波方法［J］.武漢大學學報（工學版），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 文獻標識碼：B 文章編號：1672-9900（2017）06-0078-05

2017-08-22

秦耀軍（1964-），男（漢族），山東濟寧人，高級工程師，主要從事水文地質、工程地質、環境地質、地災防治研究工作，（Tel）18561188280。

王艷肖）