本發(fā)明涉及地下水?dāng)?shù)值模擬多模型研究領(lǐng)域，具體涉及一種地下水?dāng)?shù)值模擬的多模型構(gòu)建方法。

背景技術(shù)：

地下水模擬技術(shù)已成為研究地下水環(huán)境各類問題的主要方法。地下水模擬過程首先是對研究對象的地質(zhì)和水文地質(zhì)條件高度概化并得到一個理想的物理模型，再用簡潔的數(shù)學(xué)語言去刻畫它的數(shù)量關(guān)系和空間形式，使它能夠反映研究對象的地質(zhì)、水文地質(zhì)和地下水動力特征，從而達到數(shù)值再現(xiàn)實際水流系統(tǒng)基本情況的目的。

水文地質(zhì)參數(shù)的空間變化是水流及溶質(zhì)運移模擬不確定性的主要原因。非均質(zhì)性是介質(zhì)滲透性參數(shù)的普遍特征，研究滲透性參數(shù)的非均質(zhì)特征以及由此引起的空間變異性是研究地下水滲流和溶質(zhì)運移的基礎(chǔ)。

水文地質(zhì)模擬過程中第一步并是地下水系統(tǒng)的概化，也是在模擬中最為關(guān)鍵的一步。如美國科學(xué)院院士anderson就強調(diào)“模型預(yù)報未來的失敗不是由于模型中數(shù)值或理論上的缺陷，或者確切地說，預(yù)報中的錯誤是由于概念模型有錯。王浩院士在總結(jié)地下水?dāng)?shù)值模擬的方法論時指出，地下水?dāng)?shù)值模型僅是定量刻畫地下水系統(tǒng)的工具，模擬結(jié)果的合理性和可靠性最主要取決于水文地質(zhì)概念模型的合理概化，如含水層結(jié)構(gòu)，邊界條件等，而非僅僅依靠模型方法本身。對地下水系統(tǒng)的概化從系統(tǒng)論的觀點相當(dāng)于構(gòu)建一種系統(tǒng)結(jié)構(gòu)，輸入經(jīng)過系統(tǒng)的變換而產(chǎn)生輸出，這種變化就取決與系統(tǒng)的結(jié)構(gòu)。研究表明，在同等降水條件下，不同的地下水系統(tǒng)，由于其巖層、構(gòu)造、地貌及邊界條件不同，地下水滲流場中的流向及水力梯度的變化各不相同。所以地下水模型的建立，應(yīng)該充分考慮模型概化的不確定性。

傳統(tǒng)多模型分析遵循下面的主要步驟：①考慮構(gòu)建模擬區(qū)多個可能的模型；②在相同觀測數(shù)據(jù)條件下校正這些模型；③使用某種準(zhǔn)則對模型進行排序；④去掉可能性小的模型；⑤對余下模型得到的預(yù)測值與統(tǒng)計量進行權(quán)重分析。這種多模型分析存在數(shù)據(jù)處理量大、處理步驟多的問題；且這種多模型分析對地下水模擬的高度概化與水文地質(zhì)條件及問題本身復(fù)雜性的不協(xié)調(diào)，導(dǎo)致地下水模擬不確定性問題的出現(xiàn)，引起模擬預(yù)測結(jié)果與實際情況的偏差，致使地下水環(huán)境問題的評估及治理決策存在風(fēng)險。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的是：針對現(xiàn)有技術(shù)的不足，提供一種地下水?dāng)?shù)值模擬的多模型構(gòu)建方法，以期提高地下水模擬模型精度，為地下水環(huán)境自凈模擬研究提供基礎(chǔ)模型。

為了解決上述技術(shù)問題，本發(fā)明采用的技術(shù)方案概述如下：

一種地下水?dāng)?shù)值模擬的多模型構(gòu)建方法，包括如下步驟：

(1)根據(jù)原始地形資料及現(xiàn)實測地形數(shù)據(jù)，構(gòu)建模型1#和模型2#，模型1#為考慮堆填體對原始地形的改變，模型2#為原始地形；

(2)根據(jù)收集及原位水文地質(zhì)試驗獲取的滲透系數(shù)數(shù)據(jù)，生成滲透系數(shù)對數(shù)的初始平均值及初始方差，并按am-mcmc算法生成每組樣本，每組樣本中包含對應(yīng)網(wǎng)格數(shù)的滲透系數(shù)對數(shù)值，將每組對數(shù)值進行轉(zhuǎn)化并輸入模型1#進行模擬計算，根據(jù)模擬計算的接受條件篩選含一定數(shù)據(jù)的100組參數(shù)并獲取相應(yīng)輸出數(shù)據(jù)；再將篩選好的100組參數(shù)輸入模型2#獲取相應(yīng)輸出數(shù)據(jù)；

(3)依據(jù)aicc準(zhǔn)則進行分析，基于aicc準(zhǔn)則的參數(shù)識別是利用步驟(2)得到的模型1#和模型2#的輸出數(shù)據(jù)來計算平均模型預(yù)測值及模型殘差，并通過排列模型，計算模型概率或權(quán)重來分析模型的最優(yōu)取值區(qū)間。

進一步地，所述am-mcmc算法是將參數(shù)組看成多維的向量抽取新的參數(shù)向量樣本時，第i步參數(shù)的分布為均值θi、協(xié)方差為ci的多元正態(tài)分布，其中,協(xié)方差公式為式1，在i0初始迭代中，協(xié)方差矩陣ci取固定值c0，之后自適應(yīng)更新：

式中：i0—初始化階段的樣本數(shù)；c0—初始協(xié)方差；ci—第i步參數(shù)分布的協(xié)方差；其中，ε—較小的數(shù)，確保ci不成為奇異矩陣；sd—比例因子，依賴于參數(shù)空間維度d，以確保接受率在一個合適的范圍內(nèi)；id—d維的單位矩陣，本次研究中ε＝10^-5，sd＝2.4²/d，d為參數(shù)個數(shù)；

第i+1次迭代，由公式1推出協(xié)方差公式2：

式中：ci+1—i+1步參數(shù)分布的協(xié)方差；—前i-1和i次參數(shù)的均值；a-m算法采樣具體過程如下：

①按先驗分布隨機產(chǎn)生初始樣本θ0；

②利用公式計算ci；

③產(chǎn)生參數(shù)值θ*～n(θi，ci)；

④計算接受樣本條件如式3：

式中：n—觀測數(shù)據(jù)個數(shù)；f(θi+1)—模型參數(shù)取值為θi+1得到的模擬值；y為觀測值向量；k—樣本輸出數(shù)據(jù)殘差平方均值接受范圍，取水位年紀(jì)動態(tài)變化值；

⑤重復(fù)步驟②～④，直到取得足夠多的樣本為止。

進一步地，所述aicc信息量的計算式如式4-6：

式中：si及—殘差平方和及最優(yōu)目標(biāo)函數(shù)值的殘差平方和；σi²及—殘差平方均值及最優(yōu)目標(biāo)函數(shù)值的殘差平方均值；k—待估參數(shù)個數(shù)；

得到aicc之后，用模型的aicc值減去所有備選模型中的aiccmin，計算每個模型的delta值△i，如式7：

δi＝aicci-aiccmin(式7)

最后根據(jù)delta值計算模型的后驗概率ωi，r是參加多模型分析得備選模型總數(shù)，如式8：

相對于現(xiàn)有技術(shù)，本發(fā)明所產(chǎn)生的有益效果：

1、本發(fā)明將am-mcmc算法及aicc信息量準(zhǔn)則結(jié)合構(gòu)建的多模型分析方法可綜合分析模型不確定性及參數(shù)不確定對模型精度的影響，充分考慮了模型概化的不確定性，提高了地下水模擬模型精度，為地下水環(huán)境自凈模擬研究提供了一個基礎(chǔ)模型，同時可分析各因子對模型影響的敏感程度；

2、本發(fā)明中的a-m算法的接受條件由似然函數(shù)求解的后驗概率，修正為更直接的殘差平方是否在預(yù)設(shè)的值域范圍，將傳統(tǒng)多模型分析的5個步驟縮減至3個步驟，減少了數(shù)據(jù)量，提高了模型校驗效率。

附圖說明

下面結(jié)合附圖對本發(fā)明作進一步說明。

圖1是實施例中參數(shù)在采樣過程中均值的迭代跡線圖；

圖2是實施例中參數(shù)在采樣過程中方差的迭代跡線圖；

圖3是實施例中45000個參數(shù)樣本的采樣過程中隨機樣本的參數(shù)對數(shù)值分布圖；

圖4是實施例中45000個參數(shù)樣本的采樣過程中相應(yīng)滲透系數(shù)對數(shù)值的相對頻率分布圖；

圖5是實施例中模型1#計算模型地下水位眾數(shù)、95％置信區(qū)間(陰影區(qū)域)與觀測值分布圖；

圖6是實施例中不同精度模型比例圖。

具體實施方式

下面結(jié)合附圖和具體實施方式對本發(fā)明作進一步詳細(xì)說明。本發(fā)明的實施方式包括但不限于下列實施例。

實施例

對位于成都平原區(qū)金馬河的一級階地進行地下水?dāng)?shù)值模擬的多模型構(gòu)建，包括如下步驟：

(1)模型概化及邊界條件設(shè)置

研究區(qū)位于成都平原區(qū)金馬河的一級階地，下伏含水層為第四系全新統(tǒng)沖洪積層砂卵礫石孔隙潛水(q4^al+pl)，根據(jù)鉆探成果及模擬區(qū)原位水文地質(zhì)試驗，研究區(qū)含水層厚度約為20m，東側(cè)為當(dāng)?shù)刈畹颓治g基準(zhǔn)面金馬河，地下水由北西向南東徑流。區(qū)內(nèi)1995年運營生活垃圾填埋場，2010年停止堆填，2013年采取封場措施，至今原始地形地貌已發(fā)生改變，根據(jù)收集的原始地形資料及現(xiàn)堆填區(qū)實測地形數(shù)據(jù)對比，堆填區(qū)經(jīng)過削坡、整形等封場措施后高于原始地形約9m。地表高程及坡降等因素的改變將對地下水補給條件產(chǎn)生影響，文中將以地形地貌的差異構(gòu)建模型1#和模型2#，模型基本參數(shù)取值如表1，構(gòu)建模型范圍x×y：1000×900(x為南北向，y為東西向)，網(wǎng)格為10×10m，根據(jù)水文地質(zhì)條件及鉆探成果，含水層厚度平均厚度為20m，研究區(qū)多年平均降雨量為1243mm，研究區(qū)位于岷江片區(qū)第四系沖洪積扇頂，根據(jù)水文地質(zhì)特征及多年水文資料統(tǒng)計，平均降雨入滲系數(shù)為0.186，模型設(shè)置recharge為231.27mm。東側(cè)邊界為當(dāng)?shù)刈畹颓治g基準(zhǔn)面設(shè)置為河流邊界，西側(cè)及北側(cè)根據(jù)鉆孔數(shù)據(jù)設(shè)定為定水頭邊界，計算區(qū)內(nèi)設(shè)置14個水位觀測孔。

表1模型基本參數(shù)取值

(2)am-mcmc采樣

含水介質(zhì)的非均質(zhì)問題的研究是隨著地下水隨機模擬理論提出的。在傳統(tǒng)的多孔介質(zhì)流體力學(xué)中，介質(zhì)的非均質(zhì)性決定了表征其參數(shù)的隨機性和結(jié)構(gòu)性。隨機性指參數(shù)在特定點上取值不確定；結(jié)構(gòu)性則是指同類介質(zhì)的屬性表征參數(shù)總體取值服從一定概率分布特征；結(jié)構(gòu)性是參數(shù)非均質(zhì)研究的關(guān)鍵。law首先在研究油田勘探的巖芯資料時，發(fā)現(xiàn)滲透系數(shù)服從對數(shù)正態(tài)分布，在隨后的一系列研究中^[94]均相繼證實與law相同的結(jié)論。

本實施例將表征介質(zhì)滲透性的滲透系數(shù)視為隨機變量，應(yīng)用am-mcmc算法進行采樣。將收集的模擬區(qū)周邊同類含水介質(zhì)的水文地質(zhì)試驗數(shù)據(jù)，含129個鉆孔，測得如表2所示的174組抽水試驗數(shù)據(jù)；及模擬區(qū)原位水文地質(zhì)試驗數(shù)據(jù)，含5個鉆孔，測得如表3所示的數(shù)據(jù)為先驗信息。經(jīng)統(tǒng)計，先驗信息對數(shù)取值范圍1.369～5.583，其服從正態(tài)分布。根據(jù)數(shù)據(jù)估算，參數(shù)初始樣本服從均值3.407，協(xié)方差c0為0.713的正態(tài)分布。通過a-m算法采集參數(shù)樣本，結(jié)合地下水?dāng)?shù)值模擬軟件modflow輸出模擬結(jié)果，統(tǒng)計參數(shù)樣本及模擬結(jié)果進行不確定性分析。

表2研究區(qū)同類含水介質(zhì)水文地質(zhì)試驗數(shù)據(jù)統(tǒng)計

表3研究區(qū)原位水文地質(zhì)試驗數(shù)據(jù)統(tǒng)計

本實施例中的a-m算法的接受條件由似然函數(shù)求解的后驗概率，修正為更直接的殘差平方是否在預(yù)設(shè)的值域范圍?；诖私邮軛l件的修改，首先需檢驗條件改變后，參數(shù)取值的遍歷性及收斂性。

(1)參數(shù)取值歷遍性分析

采用a-m算法對參數(shù)樣本進行采樣，根據(jù)接受條件共篩選100組參數(shù)，每組含2250個樣本值。圖1、圖2為參數(shù)在采樣過程中均值與方差的迭代跡線。當(dāng)取樣20組(樣本個數(shù)達到45000以上時)，參數(shù)的均值和方差趨于平穩(wěn)。圖3、圖3分別給出45000個參數(shù)樣本的采樣過程，樣本值遍歷參數(shù)的可能取值范圍，通過方差的自適應(yīng)更新，樣本值取樣波動逐步減弱，采樣過程基本穩(wěn)定。綜合考慮均值、方差迭代跡線和樣本采樣過程，調(diào)整接受條件后的a-m采樣方法并沒有對參數(shù)后驗分布的遍歷性及收斂性產(chǎn)生影響，更直接的接受條件可優(yōu)化樣本取樣空間，提高樣本采集效率。

(2)多模型識別分析

多模型的構(gòu)建主要目的在于分析參數(shù)不確定性與模型不確定性對模擬結(jié)果的影響程度。在剔除不符合接受條件的模型后，2組不同的概念模型分別選取100個計算模型。并依據(jù)aicc準(zhǔn)則多模型分析方法計算得到各模型的aicc值、delta值及模型的后驗概率(如表4、表5所示)。

表4模型1#排名結(jié)果(部分計算模型)

表5模型2#排名結(jié)果(部分計算模型)

aicc準(zhǔn)則分析參數(shù)不確定性對模擬預(yù)測結(jié)果的影響時，burnham和anderson建議如果最佳的模型的后驗概率超過0.9時，可視其為最佳模型來用于預(yù)測。而在地下水模擬不確定性研究中，是不易輸出如此高概率模型而獲取單個的最優(yōu)解，輸出的是一系列擬合程度相近的模型，即分析所得的是參數(shù)樣本的最優(yōu)取值區(qū)間而非唯一解。分析過程中可提供預(yù)測平均值或預(yù)測值置信區(qū)間等方法來反映預(yù)測的范圍及與參數(shù)的不確定性關(guān)系。aicc分析認(rèn)為delta值＜2的模型是較好的模型，4＜delta值＜7的模型為經(jīng)驗推薦模型，而delta值＞10的模型可以舍去。甄選模型1#中delta值小于10的計算模型，對輸出觀測孔地下水位序列值進行統(tǒng)計分析，可計算各觀測孔的眾數(shù)和置信區(qū)間。觀測孔水位眾數(shù)、95％置信區(qū)間和實際觀測值的序列如圖5所示，再進一步剔除delta值大于10的模型后，剩余計算模型的水位眾數(shù)與觀測值幾乎重合。

根據(jù)模型1#與模型2#的輸出結(jié)果排序，排名前三的參數(shù)樣本相同，排名前10的計算模型中有6組相同，表明在本文設(shè)定的不確定性條件下，較之地形條件差異的模型，參數(shù)的不確定性的敏感性更高，對模型輸出結(jié)果更具控制。雖地下水流場模擬過程中不僅會出現(xiàn)的“異參同效”，甚至存在“異參、異構(gòu)同效”，但仍能從以下兩個方面分析考慮地形變化的模型1#優(yōu)于模型2#，反映概念模型的不確定性仍對模擬結(jié)果有著明顯的影響。根據(jù)輸出結(jié)果統(tǒng)計分析，①高精度模型的出現(xiàn)頻次的不同，方差值＜1的計算模型中，模型1#與模型2#分別為6％，2％；方差值＜2的模型，模型1#與模型2#分別為65％，46％(如圖6所示)；②最優(yōu)解區(qū)間取值范圍及概率的不同，剔除delta值＞10的計算模型后，模型1#中僅保留前10個模型，累計后驗概率0.996；2#模型保留前21個模型，而其前10個模型的累計后驗概率僅為0.884(如表4、表5所示)。

如上所述即為本發(fā)明的實施例。本發(fā)明不局限于上述實施方式，任何人應(yīng)該得知在本發(fā)明的啟示下做出的結(jié)構(gòu)變化，凡是與本發(fā)明具有相同或相近的技術(shù)方案，均落入本發(fā)明的保護范圍之內(nèi)。