本發(fā)明涉及一種地下水模型評價方法，具體涉及一種基于AM嵌套抽樣算法的地下水模型評價方法。

背景技術(shù)：

近年來，數(shù)值模擬技術(shù)已成為地下水研究領(lǐng)域中一種不可或缺的方法，對于水資源評價、開發(fā)、管理與保護(hù)、地下水污染防治等問題具有重要意義。地下水模型的使用不僅可為決策者提供參考依據(jù)，也可對未來進(jìn)行預(yù)測和估計。建模方法和工具有很多，基于不同的原理或假設(shè)條件，可以建立多個不同的地下水模型。然而，選擇不同的模型對預(yù)測結(jié)果的準(zhǔn)確性有較大的影響，因此，如何評價和選擇地下水模型是當(dāng)前需要解決的問題。

邊緣似然值(綜合似然值、貝葉斯證據(jù))是評價模型表現(xiàn)或計算模型權(quán)重的重要依據(jù)，然而模型的邊緣似然值是其似然函數(shù)在復(fù)雜空間內(nèi)的高維積分，直接計算十分困難。計算邊緣似然值有多種方法，主要有：①拉普拉斯近似方法；②算數(shù)平均方法(AME)；③調(diào)和平均方法(HME)等。上述常用方法主要存在以下問題：①拉普拉斯近似方法依賴于邊緣似然值的解析形式，不適用于解析形式不存在的情況；②AME在參數(shù)先驗(yàn)分布空間內(nèi)隨機(jī)抽樣，收斂速度慢且得到的邊緣概率值偏?。虎跦ME在參數(shù)后驗(yàn)分布空間內(nèi)隨機(jī)抽樣，計算穩(wěn)定性差且容易高估邊緣似然值。

John Skilling(2006)提出了一種計算邊緣似然值的新方法：嵌套抽樣算法(Nested Sampling Algorithm，NSE)。該方法基于貝葉斯理論，其核心是將復(fù)雜的高維積分轉(zhuǎn)化為便于數(shù)值計算的一維積分。不同于AME或HME僅在先驗(yàn)或后驗(yàn)概率空間內(nèi)抽樣，也不是簡單地將先驗(yàn)與后驗(yàn)空間混合，嵌套抽樣法在抽樣的過程中由先驗(yàn)空間逐步過渡到后驗(yàn)空間，從而有效降低了從單一分布抽樣引起的邊緣似然值估計誤差。嵌套抽樣算法可以看作一種全局優(yōu)化算法，因?yàn)槠淅玫挠行?shù)集遍歷了整個先驗(yàn)分布及后驗(yàn)分布。目前，嵌套抽樣方法已經(jīng)在多個領(lǐng)域得到推廣應(yīng)用，如Elsheikh(2013)等將基于Metropolis-Hasting的嵌套抽樣算法(NSE-MH)應(yīng)用于地下水流模型的評價與不確定性分析，驗(yàn)證了嵌套抽樣算法的有效性；Liu(2016)等對NSE-MH中Metropolis-Hasting算法進(jìn)行改進(jìn)，分別應(yīng)用于線性、非線性函數(shù)的邊緣似然值計算，并與算術(shù)平均、調(diào)和平均及熱力學(xué)積分(TIE)方法的計算結(jié)果對比，驗(yàn)證了改進(jìn)后的嵌套抽樣算法的計算精度與效率。

嵌套抽樣算法(NSE)將復(fù)雜的高維積分問題轉(zhuǎn)化為一維積分問題。在具體實(shí)施過程中，嵌套抽樣算法可以分為嵌套抽樣主算法和局部限制抽樣子算法(Constrained local sampling)兩部分，主算法通過有效集迭代更新的方式實(shí)現(xiàn)嵌套抽樣，局部限制抽樣負(fù)責(zé)生成每次迭代過程所需的似然值～先驗(yàn)分布累積(L～X)樣本。局部限制抽樣通?；诟怕食闃臃椒ǎ鏜etropolis-Hasting(MH)算法等。

對于常規(guī)的基于Metropolis-Hasting算法的嵌套抽樣算法(NSE-MH)，該算法原理簡單，容易操作，但在應(yīng)用過程中存在以下問題：①NSE-MH算法的計算效率低，所需的計算量大；②NSE-MH算法的收斂速度慢，在抽樣后期MH算法需要多次迭代才能生成滿足約束條件(L_i+l＞L_i)的樣本；③NSE-MH算法在參數(shù)后驗(yàn)分布空間內(nèi)隨機(jī)抽樣，計算穩(wěn)定性穩(wěn)定性較差。因此，上述問題的存在限制了嵌套抽樣算法在模型評價中應(yīng)用和推廣。

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

發(fā)明目的：本發(fā)明的目的在于針對現(xiàn)有技術(shù)的不足，提供一種精確、高效的基于AM嵌套抽樣算法的地下水模型評價方法。

技術(shù)方案：本發(fā)明提供了一種基于AM嵌套抽樣算法的地下水模型評價方法，包括以下步驟：

(1)根據(jù)研究區(qū)的水文地質(zhì)條件，建立一組可行的概念模型M^k(k＝1，2，…，K)來表示實(shí)際地下水系統(tǒng)，K表示概念模型的數(shù)量，這些概念模型具有不同的結(jié)構(gòu)；

(2)根據(jù)研究問題選擇一組水文地質(zhì)參數(shù)作為參數(shù)向量θ并根據(jù)相關(guān)監(jiān)測資料確定其先驗(yàn)概率分布p(θ|M^k)，先驗(yàn)概率分布通常為均勻分布；

(3)從先驗(yàn)分布p(θ|M^k)中隨機(jī)生成參數(shù)向量θ的集合S＝{θ₁，θ₂，…，θ_N}作為有效集，并計算有效集中每個參數(shù)向量的聯(lián)合似然函數(shù)L(θ|D，M^k)，；

(4)確定嵌套抽樣主算法的迭代次數(shù)R，在每次迭代過程中選出有效集S中最差的參數(shù)向量作為樣本，并根據(jù)梯形公式計算邊緣似然值的增量ΔZ；

(5)在每次迭代過程中，通過基于AM算法的局部限制抽樣從先驗(yàn)分布p中生成新的參數(shù)向量θ_new作為候選樣本，以替代有效集中最差的樣本；

(6)完成R次迭代后，根據(jù)有效集S和邊緣似然值的增量ΔZ，計算各個概念模型的邊緣似然值Z；

(7)根據(jù)各個模型M^k(k＝1，2，…，K)的邊緣似然值，對各個概念模型進(jìn)行評價。

進(jìn)一步，步驟(3)計算聯(lián)合似然函數(shù)L(θ|D，M^k)：

式中，C為協(xié)方差矩陣，為單位矩陣I_d，μ為研究區(qū)地下水實(shí)測數(shù)據(jù)，Y為根據(jù)參數(shù)向量θ和模型通過數(shù)值模擬得到的數(shù)據(jù)，μ和Y是與地下水模型相關(guān)的狀態(tài)變量，如地下水水位、地下水中污染物的濃度、溫度等，n為實(shí)測值和模擬值的個數(shù)。

進(jìn)一步，步驟(4)對于第i(i＝1，…，R)次迭代，計算有效集S中最小的參數(shù)向量θ_worst及其對應(yīng)的似然函數(shù)L_worst，令L_i＝L_worst，計算先驗(yàn)分布累積X_i(X_i與有效集中參數(shù)向量的個數(shù)N和迭代次數(shù)i有關(guān))、每一次迭代中的邊緣似然值Zi以及邊緣似然值的增量ΔZ，其中Z₀＝0，L₀＝0：

進(jìn)一步，步驟(5)通過局部限制抽樣從參數(shù)先驗(yàn)分布中生成新參數(shù)向量θ_new，若L(θ_new|D，M)＞L_worst，則用θ_new取代原有θ_worst；否則，繼續(xù)從局部限制抽樣算法中生成θ_new，直至滿足L(θ_new|D，M)＞L_worst或達(dá)到人為定義的抽樣次數(shù)上限為止。

進(jìn)一步，步驟(5)基于AM算法的局部限制抽樣包括以下步驟：

①從有效集S中隨機(jī)選擇某一參數(shù)向量θ作為初始參數(shù)向量

②確定AM算法的循環(huán)次數(shù)H，對于第j(j＝1，…，H)次循環(huán)，從正態(tài)分布N(C_j)中生成新樣本ξ，計算對應(yīng)的聯(lián)合似然函數(shù)值L_ξ，其中C_j為協(xié)方差矩陣；

在T₀次迭代前取固定值C₀，之后自適應(yīng)更新協(xié)方差矩陣計算公式如下：

式中，為已有的所有參數(shù)向量的協(xié)方差矩陣；

為方便計算，可以通過遞歸公式計算C_j+1：

式中，s_d＝(2.4)²/d，d是參數(shù)的維度，堤一個大于0的常數(shù)，I_d是d維單位矩陣，和分別表示前j-1次和j次的抽樣的均值；

③若L_ξ＞L_worst，則計算接受概率否則α＝0；

④從均勻分布U(0，1)中生成隨機(jī)數(shù)u，比較u與α的大??；若u≤α則接受否則

⑤重復(fù)步驟②-④，直至生成長度為H的馬爾可夫鏈為止；令

進(jìn)一步，步驟(6)分別計算當(dāng)前有效集S中的N個參數(shù)向量θ₁，θ₂，…，θ_N對應(yīng)的似然函數(shù)L₁，L₂，…，L_N，計算得到邊緣似然值z：

有益效果：本發(fā)明將嵌套抽樣算法中的局部限制抽樣算法改進(jìn)為AM算法，將模型的邊緣似然值及后驗(yàn)概率(權(quán)重)作為評價地下水模型表現(xiàn)的指標(biāo)，根據(jù)貝葉斯分析理論及嵌套抽樣算法，將復(fù)雜且不易直接求解的高維積分邊緣似然值轉(zhuǎn)化為易于計算的一維積分，在計算地下水模型邊緣似然值的案例分析中，本方法通過AM的自適應(yīng)更新，保證了抽樣的質(zhì)量與精度，與原有的NSE-MH算法相比，在計算結(jié)果的計算效率和收斂速度方面有所提高，同時也提高了計算結(jié)果的準(zhǔn)確性和穩(wěn)定性。

附圖說明

圖1為基于AM的嵌套抽樣主算法流程圖；

圖2為研究區(qū)平面示意圖；

圖3為實(shí)施例中四個模型M1、M2、M3和M4分別采用三種方法進(jìn)行地下水流模型抽樣過程圖。

具體實(shí)施方式

下面對本發(fā)明技術(shù)方案進(jìn)行詳細(xì)說明，但是本發(fā)明的保護(hù)范圍不局限于所述實(shí)施例。

實(shí)施例：如圖1所示，一種基于AM嵌套抽樣算法的地下水模型評價方法，具體操作如下：

本實(shí)施例以三維穩(wěn)定流地下水流模型為例，研究區(qū)如圖2所示，為一矩形含水層，東西方向長為5000m，南北方向?qū)挒?000m。含水層的總厚度為60m，從上至下依次為潛水含水層、弱透水層和承壓含水層，厚度依次為35m、5m和20m。滲透系數(shù)K具有非均質(zhì)性，滲透系數(shù)隨機(jī)場用各向同性的指數(shù)方差模型描述，相關(guān)長度為200m，滲透系數(shù)K的對數(shù)(logK)的方差為1.0。各層滲透系數(shù)的平均值依次為1.0m/d、0.1m/d和5.0m/d。此外，含水層水平方向的滲透系數(shù)是垂向滲透系數(shù)的10倍。

研究區(qū)東側(cè)為一河流邊界，僅切穿潛水含水層，河水位為35m，河床底板高程為30m，河床水力傳導(dǎo)系數(shù)為20m²/d；西側(cè)為一定水頭邊界，水位為56m；右側(cè)有一排水渠，渠道底板高程為45m，渠的水力傳導(dǎo)系數(shù)為20m²/d。此外，研究區(qū)的南部、北部和底部均為隔水邊界。潛水含水層接受大氣降水的均勻補(bǔ)給，降水量為9.0×10^-4m/d，入滲補(bǔ)給系數(shù)為0.15。研究區(qū)內(nèi)共有5口抽水井，總抽水量為1250m³/d(每口井的抽水量均為250m³/d)；區(qū)內(nèi)有32口地下水位監(jiān)測井(第1層和第3層中各16口)。

考慮到模型結(jié)構(gòu)的不確定性，本次研究采用四個結(jié)構(gòu)不同的概念模型(M1、M2、M3和M4)作為對未知地下水系統(tǒng)的近似。該模型結(jié)構(gòu)的不確定性主要體現(xiàn)為對中間弱透水層的刻畫，因?yàn)樵趯?shí)際水文地質(zhì)調(diào)查中往往很難查明弱透水層的位置及厚度，而弱透水層又對地下水系統(tǒng)的模擬結(jié)果有很大的影響。4個概念模型的空間大小與真實(shí)情況完全相同(5000m×3000m×60m)。模型M1假設(shè)只有一層潛水含水層；模型M2假設(shè)存在潛水含水層和承壓含水層，厚度分別為35m和25m；模型M3、M4的結(jié)構(gòu)與真實(shí)情況相同，M3中含水層的厚度分別為35m、3m和22m；M4含水層的厚度分別為35m、7m和18m。

本次研究5個水文地質(zhì)參數(shù)，分別為第一層的入滲補(bǔ)給系數(shù)、定水頭邊界水頭、河床水力傳導(dǎo)系數(shù)、第一層的滲透系數(shù)隨機(jī)場的方差和相關(guān)長度。其他模型邊界條件與真實(shí)模型相同。5個參數(shù)的先驗(yàn)信息均為均勻分布，其取值范圍見表1：

表1地下水流模型中參數(shù)的取值范圍

定義有效集中樣本的個數(shù)N＝25，根據(jù)公式(1)計算每個樣本的邊緣似然值L，其中n＝32，μ為地下水監(jiān)測井中的水位實(shí)測值，Y是將樣本代入地下水流模型中得到的模擬值。本例中水流模型基于MODFLOW-2005建立并求解。

在本例中，確定嵌套抽樣主算法的迭代次數(shù)R＝250，在第i(i＝1，…，250)次迭代過程中，選出有效集中的最小的參數(shù)向量θ_worst及其對應(yīng)的似然函數(shù)L_worst，令L_i＝L_worst，根據(jù)公式(2)計算先驗(yàn)分布累積X_i，根據(jù)梯形公式(3)計算邊緣似然值的增量ΔZ和每一次迭代中的邊緣似然值Z_i。

在第i(i＝1，…，250)次迭代中，通過局部限制抽樣(AM算法)從參數(shù)先驗(yàn)分布中生成新參數(shù)向量θ_new，若L(θ_new|D，M)＞L_worst，則用θ_new取代原有θ_worst；否則，繼續(xù)從局部限制抽樣算法中生成θ_new，直至滿足L(θ_new|D，M)＞L_worst或達(dá)到人為定義的抽樣次數(shù)上限為止，本例中設(shè)置的抽樣次數(shù)上限為200。

在基于AM算法的局部限制抽樣中，從有效集中隨機(jī)選擇某一參數(shù)向量θ作為初始參數(shù)向量確定AM算法的循環(huán)次數(shù)H＝100，對于第j(j＝1，…，100)次循環(huán)，從正態(tài)分布N(C_j)中生成新樣本ξ，根據(jù)公式(1)計算對應(yīng)的聯(lián)合似然函數(shù)值L_ξ，其中C_j為協(xié)方差矩陣，在T₀次迭代前取固定值C₀，本例中T₀＝20，C₀為單位矩陣。隨后根據(jù)公式(4)和(5)自適應(yīng)更新協(xié)方差矩陣；若L_ξ＞L_worst，則計算接受概率α，否則α＝0；從均勻分布U(0，1)中生成隨機(jī)數(shù)u，比較u與α的大小；若u≤α則接受否則直至生成長度為100的馬爾可夫鏈為止；令

完成250次迭代后，根據(jù)有效集和邊緣似然值的增量ΔZ，計算各個概念模型的邊緣似然值Z；分別計算當(dāng)前有效集中的25個參數(shù)向量θ₁，θ₂，…，θ₂₅對應(yīng)的似然函數(shù)L₁，L₂，…，L₂₅，根據(jù)公式(6)計算得到邊緣似然值Z。

為驗(yàn)證NSE-AM算法的有效性，同時采用算術(shù)平均方法(AME)和NSE-MH算法求解。本次案例分析中算術(shù)平均方法的樣本數(shù)量為50萬，分別計算4個模型的邊緣似然值，作為對應(yīng)的參考值；NSE-MH的參數(shù)取值為：N＝25，R＝250，局部限制抽樣算法中生成的馬爾可夫鏈長為100。使用NSE-AM和NSE-MH對4個模型各進(jìn)行似然值估計10次。由于嵌套抽樣算法的各次計算結(jié)果不盡相同，本次研究中取10次計算結(jié)果的平均值作為最終結(jié)果，從而評價嵌套抽樣算的綜合表現(xiàn)。根據(jù)三種方法的迭代過程繪制邊緣似然值隨抽樣次數(shù)的變化趨勢圖，如圖3所示；三種方法收斂后計算得到的邊緣似然值如表2所示：

表2兩種方法計算水流模型的邊緣似然值

從中可以得到以下結(jié)論：(1)NSE-AM算法的結(jié)果準(zhǔn)確、計算效率高，AME方法通常需要數(shù)十萬次的樣本數(shù)量(即地下水流模型的運(yùn)行次數(shù))才能達(dá)到穩(wěn)定，需要較長的計算時間，NSE-MH算法也需要5萬次以上的模型運(yùn)行才能收斂，而NSE-AM算法經(jīng)過大約2萬次的模型運(yùn)行即收斂，且最終得到的結(jié)果與AME方法較為接近；(2)NSE-AM算法的適應(yīng)性好，對于不同的概念模型都能在有限的模型運(yùn)行次數(shù)中達(dá)到穩(wěn)定，而AME方法對于個別與真實(shí)情況相差較大的模型，如M1經(jīng)過50萬次模型運(yùn)行也未達(dá)到穩(wěn)定，而NSE-AM算法在2萬次模型運(yùn)行后基本達(dá)到穩(wěn)定；(3)根據(jù)4個模型計算得到的邊緣似然值，從小到大的順序依次為：M1＜M2＜M4＜M3，說明模型M3和M4明顯優(yōu)于模型M1和M2。這表明在考慮相同參數(shù)不確定性的情況下，模型的邊緣似然值與模型結(jié)構(gòu)概化的合理程度有關(guān)，即模型的邊緣似然值越大，模型結(jié)構(gòu)概化越合理，概念模型越接近于真實(shí)模型。