王科奇, 劉佩貴, 陶月贊, 李 飛, 樊慧慧

(1.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥 230009; 2.合肥工業(yè)大學(xué) 資源與環(huán)境工程學(xué)院,安徽 合肥 230009)

礦坑涌水量預(yù)測(cè)模型參數(shù)靈敏度分析

王科奇1, 劉佩貴1, 陶月贊1, 李 飛1, 樊慧慧2

(1.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥 230009; 2.合肥工業(yè)大學(xué) 資源與環(huán)境工程學(xué)院,安徽 合肥 230009)

礦區(qū)水文地質(zhì)條件一般較復(fù)雜,參數(shù)的空間變異性強(qiáng),致使礦坑涌水量預(yù)測(cè)模型中的參數(shù)具有一定的不確定性,從而影響預(yù)測(cè)結(jié)果的可靠度。為定量確定影響預(yù)測(cè)結(jié)果的主要因子,文章以某金屬礦區(qū)的礦坑涌水量計(jì)算模型為研究對(duì)象,運(yùn)用擴(kuò)展傅里葉幅度靈敏度檢驗(yàn)法(extended Fourier amplitude sensitivity test,EFAST)對(duì)模型參數(shù)進(jìn)行了全局靈敏度分析,得到相關(guān)參數(shù)的一階靈敏度和總靈敏度,并分別與局部靈敏度分析法和Morris法的評(píng)價(jià)結(jié)果進(jìn)行了詳細(xì)對(duì)比分析。研究結(jié)果表明,采用EFAST法開(kāi)展礦坑涌水量預(yù)測(cè)模型的參數(shù)靈敏度分析是切實(shí)可行的,也是較為有效的方法之一,礦坑涌水量預(yù)測(cè)結(jié)果對(duì)滲透系數(shù)最為敏感。

礦坑涌水量;滲透系數(shù);靈敏度分析;擴(kuò)展傅里葉幅度靈敏度檢驗(yàn)法(EFAST);Morris法

礦坑涌水量的準(zhǔn)確預(yù)測(cè)是礦區(qū)水文地質(zhì)工作的重要內(nèi)容之一,其預(yù)測(cè)方法主要有解析法、數(shù)值法、水文地質(zhì)比擬法等,這些方法均要借助數(shù)學(xué)模型計(jì)算涌水量,而影響其計(jì)算精度的關(guān)鍵要素之一是參數(shù)的合理取值。然而受成礦條件等因素的影響,大多數(shù)礦區(qū)水文地質(zhì)條件較為復(fù)雜,含水層空間變異性強(qiáng),此外還受到勘察程度、認(rèn)知水平等條件的限制,礦坑涌水量預(yù)測(cè)模型中的相關(guān)參數(shù)往往存在很大的不確定性,給參數(shù)的取值造成了一定的困難,也嚴(yán)重影響了涌水量預(yù)測(cè)結(jié)果的可靠度。研究預(yù)測(cè)結(jié)果對(duì)模型中重要參數(shù)的響應(yīng)程度(即參數(shù)的靈敏度分析),是解決該問(wèn)題的有效途徑之一。

靈敏度分析方法可以定性或者定量分析模型輸入?yún)?shù)的不確定性對(duì)模型結(jié)果產(chǎn)生的影響程度,合理確定主要與次要因子。靈敏度分析可分為局部靈敏度分析(local sensitivity analysis)和全局靈敏度分析[1](global sensitivity analysis)2類(lèi)。其中,局部靈敏性分析因具有操作簡(jiǎn)便、易于理解等優(yōu)點(diǎn),在地下水?dāng)?shù)學(xué)模型中應(yīng)用較廣泛[2-4],但該方法未考慮參數(shù)間的耦合作用,其分析結(jié)果有一定的片面性。全局靈敏度分析則突破了該局限,且可以提供更加豐富的信息,分析方法主要有回歸法[5]、Morris法[6]、Sobol’法[7]、傅里葉幅度靈敏度檢驗(yàn)法[8](Fourier amplitude sensitivity test,FAST)以及擴(kuò)展傅里葉幅度靈敏度檢驗(yàn)法[9](extended FAST,EFAST)等。其中,Morris法具有計(jì)算簡(jiǎn)單、計(jì)算成本較低的優(yōu)點(diǎn),在地下水?dāng)?shù)學(xué)模型分析方面應(yīng)用較多[10-12],但該方法是一種定性全局靈敏度分析方法,不能對(duì)參數(shù)的影響大小進(jìn)行定量化,且參數(shù)的初始值和擾動(dòng)值會(huì)對(duì)結(jié)果產(chǎn)生較大影響[11,13]。而EFAST法是一種基于方差分析的全局靈敏度分析方法,能夠定量地獲得參數(shù)的靈敏度,與其他方法相比可獲得的信息多,相對(duì)較優(yōu)。因此,本文以丘陵區(qū)某金屬礦山的礦坑涌水量預(yù)測(cè)模型為例,選取模型中的相關(guān)參數(shù),采用EFAST法分析水量預(yù)測(cè)結(jié)果對(duì)各個(gè)參數(shù)的敏感程度,以期為提高預(yù)測(cè)結(jié)果的可靠度提供科學(xué)依據(jù)。

1 傅里葉幅度靈敏度擴(kuò)展檢驗(yàn)法

EFAST法是一種基于方差分析的靈敏度分析方法,認(rèn)為參數(shù)的靈敏度可以用模型結(jié)果的方差來(lái)反映,模型結(jié)果的方差由各參數(shù)及參數(shù)間的相互作用所導(dǎo)致[14]。該法既具有FAST方法計(jì)算參數(shù)一階靈敏度的高效性,又與Sobol’法一樣有計(jì)算參數(shù)的總靈敏度的能力。

假定模型為Y=f(X),模型的輸入?yún)?shù)為X(x1,x2,…,xn),選擇適當(dāng)?shù)乃阉骱瘮?shù)Gi將模型Y=f(x1,x2,…,xn)轉(zhuǎn)化為Y=f(s),即實(shí)現(xiàn)了多維到一維的轉(zhuǎn)換。取

xi(s)=G(sinwis),i=1,2，…,n

(1)

其中,s為標(biāo)量,且-∞

通過(guò)對(duì)f(s)進(jìn)行傅里葉級(jí)數(shù)展開(kāi)可得到:

y=f(s)=∑|Ajcosjs+Bjcosjs|

(2)

其中, +∞

由參數(shù)xi的不確定性所引起的模型輸出結(jié)果方差Vi為wi的整數(shù)倍振幅平方和[15]，即

(3)

其中,Z0=Z-{0}。

模型輸出結(jié)果的總方差V為所有頻率的振幅平方和:

(4)

參數(shù)xi的一階全局靈敏度Si定義為:

Si=Vi/V

(5)

一階全局靈敏度能夠反映參數(shù)自身的不確定性對(duì)模型結(jié)果的直接影響程度。

模型輸出結(jié)果的方差V分解[16]為:

(6)

對(duì)模型結(jié)果的各高階方差進(jìn)行歸一化處理可以得到參數(shù)的各階靈敏度,例如:

Sij=Vij/V,S12…n=V12…n/V

(7)

針對(duì)一個(gè)多參數(shù)的非線性模型,參數(shù)xi的總?cè)朱`敏度為參數(shù)各階靈敏度的總和[7],即

STi=Si+Sij+Sijm+…+S12…i…n

(8)

參數(shù)的總靈敏度STi不僅包含了參數(shù)xi自身的不確定性對(duì)模型輸出結(jié)果方差的影響,同時(shí)也包含了xi與其余參數(shù)間的耦合作用對(duì)模型結(jié)果的影響。參數(shù)xi的總靈敏度STi與一階靈敏度Si的差值STi-Si可以表示該參數(shù)與其他參數(shù)之間相互作用對(duì)模型結(jié)果的間接影響;若差值較小,則說(shuō)明參數(shù)間的交互作用不明顯[17]。

2 實(shí)例研究

2.1 研究區(qū)概況

研究區(qū)為低山丘陵地形,屬亞熱帶季風(fēng)氣候區(qū),多年平均降水量1 216.2 mm,年平均蒸發(fā)量1 497.5 mm。礦區(qū)含水巖組主要有松散巖類(lèi)孔隙含水巖組,火山碎屑巖、熔巖類(lèi)空洞-裂隙含水巖組,碎屑巖、碳酸鹽巖類(lèi)巖溶-裂隙含水巖組及巖漿巖類(lèi)裂隙含水巖組。其中,碎屑巖、碳酸鹽巖類(lèi)巖溶-裂隙含水巖組為礦床的含礦層位,礦床內(nèi)各含水巖組間無(wú)明顯的隔水層存在,天然狀態(tài)下有一定的水力聯(lián)系,具有統(tǒng)一的水動(dòng)力場(chǎng)。降雨入滲為該礦床的主要充水因素,由風(fēng)化裂隙向深部補(bǔ)給。

2.2 參數(shù)靈敏度分析

礦體含水層為裂隙含水巖組,水力性質(zhì)為承壓水,并且無(wú)明顯的側(cè)向隔水介質(zhì),側(cè)向水力聯(lián)系較好;在開(kāi)采過(guò)程中,隨著排水疏干,地下水位會(huì)降到含水層頂板以下,水力性質(zhì)轉(zhuǎn)為承壓-潛水水流。因此選用均質(zhì)無(wú)限補(bǔ)給邊界條件下承壓轉(zhuǎn)無(wú)壓完整井的穩(wěn)定流涌水量計(jì)算公式,即

R0=R+r0,S=H-h,

(9)

其中,Q為礦坑涌水量預(yù)測(cè)值;K為滲透系數(shù);H為抽水前承壓含水層從底板起算的水頭值;M為承壓含水層厚度;h為抽水后承壓含水層從底板起算水頭值;S為水位降深;R0為礦井開(kāi)采區(qū)引用影響半徑;r0為礦井開(kāi)采區(qū)引用半徑;R為單井影響半徑;F為預(yù)計(jì)礦井開(kāi)采面積。

根據(jù)勘察成果及礦區(qū)開(kāi)展的相關(guān)水文地質(zhì)工作,確定參數(shù)取值范圍如下：K為0.42～0.65 m/d，M為419.56～546.43 m，H為557.86～632.39 m，r0為313.37～368.67 m。因資料有限,假定參數(shù)均服從均勻分布。

目前對(duì)于地下水?dāng)?shù)學(xué)模型,通常采用局部靈敏度法進(jìn)行靈敏度分析。本文首先采用差分法對(duì)礦坑涌水量預(yù)測(cè)模型進(jìn)行局部靈敏度分析,以各參數(shù)的平均值為基準(zhǔn)值,通過(guò)逐一改變參數(shù)的變化幅度來(lái)對(duì)模型進(jìn)行局部靈敏度分析,結(jié)果如圖1所示。

圖1 局部靈敏度系數(shù)法評(píng)價(jià)結(jié)果

由圖1可知,通過(guò)局部靈敏度分析法得到的參數(shù)靈敏度排序依次為|H|>|K|>|r0|>|M|。

本文借助Simlab軟件實(shí)現(xiàn)樣本的抽樣和結(jié)果的全局靈敏度分析,根據(jù)EFAST法的基本原理每個(gè)參數(shù)至少要采樣65次,本次試驗(yàn)設(shè)置樣本數(shù)為1 000。利用軟件的內(nèi)部函數(shù)編輯器輸入該數(shù)學(xué)模型,直接調(diào)用抽樣結(jié)果即可獲得模型的輸出結(jié)果。對(duì)模型結(jié)果進(jìn)行靈敏度分析,可得各輸入?yún)?shù)的一階靈敏度和總靈敏度,見(jiàn)表1所列。

表1 參數(shù)的一階靈敏度和總靈敏度

對(duì)比上述2種方法的分析結(jié)果可知,EFAST法得到的參數(shù)靈敏度排序與局部靈敏度分析法并不一致。這是因?yàn)榫植快`敏度分析方法是在線性模型的基礎(chǔ)上發(fā)展而來(lái)的,對(duì)于非線性模型，該方法只能分析因子對(duì)模型局部的影響。全局靈敏度分析方法能夠分析參數(shù)的變化范圍對(duì)模型結(jié)果的影響,與局部靈敏度分析方法相比,對(duì)于非線性模型其結(jié)果更趨于實(shí)際,這也正說(shuō)明在非線性模型靈敏度分析中采用全局方法的重要性。

為了與EFAST法作對(duì)照，本文采用在地下水?dāng)?shù)學(xué)模型中應(yīng)用較為廣泛的全局靈敏度分析法Morris法[10],利用Simlab軟件對(duì)該數(shù)學(xué)模型進(jìn)行靈敏度分析,獲得各參數(shù)的靈敏度分析結(jié)果,如圖2所示。圖2中,σ為各參數(shù)單位效應(yīng)的標(biāo)準(zhǔn)差;μ*為各參數(shù)單位效應(yīng)的均值。

圖2 基于Morris法的靈敏度分析結(jié)果

2.3 結(jié)果分析與討論

在局部靈敏度分析法中，某一特定參數(shù)k的靈敏度系數(shù)S是通過(guò)只改變參數(shù)k的值,使其由ak變化為ak+Δak,相應(yīng)模型結(jié)果由yi(ak)變化為yi(ak+Δak),由(10)式獲得S的近似值,即

(10)

由(10)式可知局部靈敏度分析法能夠在線性模型中獲得正確的結(jié)果,而對(duì)于非線性模型,參數(shù)在基準(zhǔn)值處變動(dòng)過(guò)大時(shí)就不能提供有效分析結(jié)果,這正是局部靈敏度分析法獲得的參數(shù)靈敏度排序與EFAST法不同的原因。

利用EFAST法能夠定量獲得各參數(shù)的一階靈敏度及總靈敏度,由表1數(shù)據(jù)可知,參數(shù)的總靈敏度與一階靈敏度之間的差值較小,說(shuō)明參數(shù)之間的耦合作用對(duì)結(jié)果的影響程度較小;滲透系數(shù)K、水頭高度H及含水層厚度M這3個(gè)參數(shù)的靈敏度相對(duì)較大,礦坑涌水量預(yù)測(cè)值對(duì)這3個(gè)參數(shù)較為敏感;其中K的不確定性對(duì)礦坑涌水量影響最大,一階靈敏度高達(dá)0.772 1。

Morris法能夠獲得各參數(shù)單位效應(yīng)的均值μ*和標(biāo)準(zhǔn)差σ,其中μ*可以表征參數(shù)靈敏度的相對(duì)大小,值越大,說(shuō)明參數(shù)的敏感性越強(qiáng);σ表征參數(shù)間的相互作用程度[18]。由圖2可得參數(shù)的靈敏度大小依次為|K|>|H|>|M|>|r0|。EFAST法獲得的一階靈敏度的排序與Morris法的排序結(jié)果是一致的,間接證明了EFAST法分析結(jié)果的可靠性,也表明應(yīng)用該方法在礦坑涌水量預(yù)測(cè)模型參數(shù)靈敏度分析方面的可行性。

3 結(jié) 論

本文基于EFAST法,以某礦區(qū)的礦坑涌水量計(jì)算模型為例,進(jìn)行了參數(shù)的靈敏度分析,分析結(jié)果表明礦區(qū)含水層的滲透系數(shù)對(duì)水量預(yù)測(cè)結(jié)果影響最大,即表明滲透系數(shù)最為敏感。為了減小礦坑涌水量預(yù)測(cè)的不確定性進(jìn)而獲得更高精度的的預(yù)測(cè)值,在勘察工作中應(yīng)盡量提高滲透系數(shù)的精度。

對(duì)于非線性模型,全局靈敏度分析結(jié)果比局部靈敏度的分析結(jié)果更趨于實(shí)際,EFAST法是定量刻畫(huà)全局靈敏度的分析方法之一,不僅克服了局部靈敏度分析方法的片面性,而且能將參數(shù)的敏感程度定量化,而不是像Morris法僅能獲得參數(shù)的靈敏性排序,提供了更綜合全面的信息。

本文選取的是解析模型,后續(xù)研究將進(jìn)行多參數(shù)數(shù)值模型的定量全局靈敏度分析。

[1] 雷曉輝,田雨,賈仰文,等.WEP 模型全局參數(shù)敏感性分析及其在漢江上游流域的應(yīng)用[J].水文,2010,30(6):14-18,41.

[2] 束龍倉(cāng),朱元生，孫慶義，等.地下水允許開(kāi)采量確定的風(fēng)險(xiǎn)分析[J].水利學(xué)報(bào),2000,31(3):77-81.

[3] 束龍倉(cāng),朱元生,孫慶義,等.地下水資源評(píng)價(jià)結(jié)果的可靠性探討 [J].水科學(xué)進(jìn)展,2000,11(1):21-24.

[4] 李森,陳家軍,葉慧海,等.地下水流數(shù)值模擬中隨機(jī)因素的靈敏度分析[J].水利學(xué)報(bào),2006,37(8):977-984.

[5] MCKAY M D,BECKMAN R J,CONOVER W J.Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J].Technometrics,1979,21(2):239-245.

[6] MORRIS M D.Factorial sampling plans for preliminary computational experiments[J].Technometrics,1991,33(2):161-174.

[7] SOBOL’ I M.Sensitivity analysis for nonlinear mathematical models[J].Mathematical Modeling and Computational Experiment,1993,1(4):407-414.

[8] CUKIER R I,FORTUIN C M,SHULER K E,et al.Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients：Ⅰ. theory[J].The Journal of Chemical Physics,1973,59(8):3873-3878.

[9] SALTELLI A,TARANTOLA S,CHAN K P S.A quantitative model-independent method for global sensitivity analysis of model output[J].Technometrics,1999,41(1):39-56.

[10] 束龍倉(cāng),王茂枚,劉瑞國(guó),等.地下水?dāng)?shù)值模擬中的參數(shù)靈敏度分析[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2007,35(5):491-495.

[11] 魯程鵬,束龍倉(cāng),劉麗紅,等.基于靈敏度分析的地下水?dāng)?shù)值模擬精度適應(yīng)性評(píng)價(jià)[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,38(1):26-30.

[12] 劉瑤林,劉國(guó)東,徐濤,等.多層含水層地下水?dāng)?shù)值模型參數(shù)靈敏度分析[J].環(huán)境科學(xué)與技術(shù),2014,37(12):33-37,44.

[13] 張麗華,蘇小四,孟祥菲,等.地下水流數(shù)值模擬參數(shù)全局靈敏度分析[J].中國(guó)農(nóng)村水利水電,2014(8):92-97.

[14] 吳錦,余福水,陳仲新,等.基于 EPIC 模型的冬小麥生長(zhǎng)模擬參數(shù)全局敏感性分析[J].農(nóng)業(yè)工程學(xué)報(bào),2009,25(7):136-142.

[15] 王建棟,郭維棟,李紅祺.拓展傅里葉幅度敏感性檢驗(yàn) (EFAST) 在陸面過(guò)程模式中參數(shù)敏感性分析的應(yīng)用探索[J].物理學(xué)報(bào),2013,62(5):35-41.

[16] SALTELLI A,BOLADO R.An alternative way to compute Fourier amplitude sensitivity test (FAST)[J].Computational Statistics & Data Analysis,1998,26(4):445-460.

[17] 余衍然,李成,姚林泉,等.基于傅里葉幅值檢驗(yàn)擴(kuò)展法的軌道車(chē)輛垂向模型全局靈敏度分析[J].振動(dòng)與沖擊,2014,33(6):77-81.

[18] 宋明丹,馮浩,李正鵬,等.基于 Morris 和 EFAST 的 CERES-Wheat 模型敏感性分析[J].農(nóng)業(yè)機(jī)械學(xué)報(bào),2014,45(10):124-131,166.

(責(zé)任編輯 張淑艷)

Sensitivity analysis for forecasting model parameters of mine water inflow

WANG Keqi1, LIU Peigui1, TAO Yuezan1, LI Fei1, FAN Huihui2

(1.School of Civil and Hydraulic Engineering, Hefei University of Technology, Hefei 230009, China; 2.School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China)

The uncertainties of mine water inflow forecasting model parameters have often been witnessed due to the complex hydrogeological conditions and high spatial variation of parameters of mining area, thus affecting the reliability of predictions. In order to quantitatively determine the main factors influencing the size of the water flow, taking the calculation model of pit water inflow of a metal mining area as the research object, the global sensitivity analysis on model parameters with the application of extended Fourier amplitude sensitivity test(EFAST) method is utilized. The first-order sensitivity of relevant parameters and the total sensitivity are obtained. In addition, the results are also compared with the parameter sensitivity sequencing obtained by the Morris method and local sensitivity analysis. The results show that the EFAST method is practical and feasible for the parameter sensitivity analysis of mine water inflow prediction model and it is one of the most effective methods, and the prediction result is most sensitive to the hydraulic conductivity.

mine water inflow; hydraulic conductivity; sensitivity analysis; extended Fourier amplitude sensitivity test(EFAST) method; Morris method

2015-09-25；

2016-02-16

國(guó)家自然科學(xué)基金資助項(xiàng)目(51309071；51509064)

王科奇(1991-),男,河南南陽(yáng)人,合肥工業(yè)大學(xué)碩士生; 陶月贊(1964-),男,安徽巢湖人,博士,合肥工業(yè)大學(xué)教授,博士生導(dǎo)師.

10.3969/j.issn.1003-5060.2017.04.014

P641.41

A

1003-5060(2017)04-0502-04