礦坑涌水量預測模型參數靈敏度分析

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

(1.合肥工業大學 土木與水利工程學院,安徽 合肥 230009; 2.合肥工業大學 資源與環境工程學院,安徽 合肥 230009)

礦坑涌水量預測模型參數靈敏度分析

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

(1.合肥工業大學 土木與水利工程學院,安徽 合肥 230009; 2.合肥工業大學 資源與環境工程學院,安徽 合肥 230009)

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

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

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

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

1 傅里葉幅度靈敏度擴展檢驗法

EFAST法是一種基于方差分析的靈敏度分析方法,認為參數的靈敏度可以用模型結果的方差來反映,模型結果的方差由各參數及參數間的相互作用所導致[14]。該法既具有FAST方法計算參數一階靈敏度的高效性,又與Sobol’法一樣有計算參數的總靈敏度的能力。

假定模型為Y=f(X),模型的輸入參數為X(x1,x2,…,xn),選擇適當的搜索函數Gi將模型Y=f(x1,x2,…,xn)轉化為Y=f(s),即實現了多維到一維的轉換。取

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

(1)

其中,s為標量,且-∞

通過對f(s)進行傅里葉級數展開可得到:

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

(2)

其中, +∞

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

(3)

其中,Z0=Z-{0}。

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

(4)

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

Si=Vi/V

(5)

一階全局靈敏度能夠反映參數自身的不確定性對模型結果的直接影響程度。

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

(6)

對模型結果的各高階方差進行歸一化處理可以得到參數的各階靈敏度,例如:

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

(7)

針對一個多參數的非線性模型,參數xi的總全局靈敏度為參數各階靈敏度的總和[7],即

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

(8)

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

2 實例研究

2.1 研究區概況

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

2.2 參數靈敏度分析

礦體含水層為裂隙含水巖組,水力性質為承壓水,并且無明顯的側向隔水介質,側向水力聯系較好;在開采過程中,隨著排水疏干,地下水位會降到含水層頂板以下,水力性質轉為承壓-潛水水流。因此選用均質無限補給邊界條件下承壓轉無壓完整井的穩定流涌水量計算公式,即

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

(9)

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

根據勘察成果及礦區開展的相關水文地質工作,確定參數取值范圍如下：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。因資料有限,假定參數均服從均勻分布。

目前對于地下水數學模型,通常采用局部靈敏度法進行靈敏度分析。本文首先采用差分法對礦坑涌水量預測模型進行局部靈敏度分析,以各參數的平均值為基準值,通過逐一改變參數的變化幅度來對模型進行局部靈敏度分析,結果如圖1所示。

圖1 局部靈敏度系數法評價結果

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

本文借助Simlab軟件實現樣本的抽樣和結果的全局靈敏度分析,根據EFAST法的基本原理每個參數至少要采樣65次,本次試驗設置樣本數為1 000。利用軟件的內部函數編輯器輸入該數學模型,直接調用抽樣結果即可獲得模型的輸出結果。對模型結果進行靈敏度分析,可得各輸入參數的一階靈敏度和總靈敏度,見表1所列。

表1 參數的一階靈敏度和總靈敏度

對比上述2種方法的分析結果可知,EFAST法得到的參數靈敏度排序與局部靈敏度分析法并不一致。這是因為局部靈敏度分析方法是在線性模型的基礎上發展而來的,對于非線性模型，該方法只能分析因子對模型局部的影響。全局靈敏度分析方法能夠分析參數的變化范圍對模型結果的影響,與局部靈敏度分析方法相比,對于非線性模型其結果更趨于實際,這也正說明在非線性模型靈敏度分析中采用全局方法的重要性。

為了與EFAST法作對照，本文采用在地下水數學模型中應用較為廣泛的全局靈敏度分析法Morris法[10],利用Simlab軟件對該數學模型進行靈敏度分析,獲得各參數的靈敏度分析結果,如圖2所示。圖2中,σ為各參數單位效應的標準差;μ*為各參數單位效應的均值。

圖2 基于Morris法的靈敏度分析結果

2.3 結果分析與討論

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

(10)

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

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

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

3 結 論

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

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

本文選取的是解析模型,后續研究將進行多參數數值模型的定量全局靈敏度分析。

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

[2] 束龍倉,朱元生，孫慶義，等.地下水允許開采量確定的風險分析[J].水利學報,2000,31(3):77-81.

[3] 束龍倉,朱元生,孫慶義,等.地下水資源評價結果的可靠性探討 [J].水科學進展,2000,11(1):21-24.

[4] 李森,陳家軍,葉慧海,等.地下水流數值模擬中隨機因素的靈敏度分析[J].水利學報,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] 束龍倉,王茂枚,劉瑞國,等.地下水數值模擬中的參數靈敏度分析[J].河海大學學報(自然科學版),2007,35(5):491-495.

[11] 魯程鵬,束龍倉,劉麗紅,等.基于靈敏度分析的地下水數值模擬精度適應性評價[J].河海大學學報(自然科學版),2010,38(1):26-30.

[12] 劉瑤林,劉國東,徐濤,等.多層含水層地下水數值模型參數靈敏度分析[J].環境科學與技術,2014,37(12):33-37,44.

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

[14] 吳錦,余福水,陳仲新,等.基于 EPIC 模型的冬小麥生長模擬參數全局敏感性分析[J].農業工程學報,2009,25(7):136-142.

[15] 王建棟,郭維棟,李紅祺.拓展傅里葉幅度敏感性檢驗 (EFAST) 在陸面過程模式中參數敏感性分析的應用探索[J].物理學報,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] 余衍然,李成,姚林泉,等.基于傅里葉幅值檢驗擴展法的軌道車輛垂向模型全局靈敏度分析[J].振動與沖擊,2014,33(6):77-81.

[18] 宋明丹,馮浩,李正鵬,等.基于 Morris 和 EFAST 的 CERES-Wheat 模型敏感性分析[J].農業機械學報,2014,45(10):124-131,166.

(責任編輯 張淑艷)

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

國家自然科學基金資助項目(51309071；51509064)

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

10.3969/j.issn.1003-5060.2017.04.014

P641.41

A

1003-5060(2017)04-0502-04