魯程鵬,束龍倉,劉麗紅,劉佩貴,董貴明

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室,江蘇南京 210098;2.合肥工業(yè)大學(xué)土木與水利工程學(xué)院,安徽合肥 230009;3.中國礦業(yè)大學(xué)資源與地球科學(xué)學(xué)院,江蘇 徐州 221116)

隨著地下水流基礎(chǔ)理論的完善和計算機(jī)技術(shù)的不斷發(fā)展,地下水?dāng)?shù)值模擬技術(shù)日趨成熟,地下水?dāng)?shù)值模型已成為地下水科學(xué)發(fā)展中一個強有力的工具[1].數(shù)值模型對定解問題的邊界條件以及初始條件的要求較解析模型低,適應(yīng)于各種復(fù)雜的地下水流系統(tǒng)問題的求解.數(shù)學(xué)理論已證明,在符合一定的時間、空間步長要求下,數(shù)值解趨近于解析解[2].但是,通過對自然條件的概化得到的數(shù)值模型的解究竟與真實情形相差多少,還是必須通過與實測資料的對比來定量評價模擬結(jié)果的優(yōu)劣.地下水領(lǐng)域使用最為廣泛的觀測資料是地下水位觀測資料.對數(shù)值模型而言,模擬結(jié)果的精度適應(yīng)性及參數(shù)靈敏度是極其重要的2個方面.模擬精度的適應(yīng)性是指模擬精度隨參數(shù)變化而做出的響應(yīng);參數(shù)靈敏度是指模型輸出在參數(shù)攝動條件下的響應(yīng).兩者存在以下差別:一是參數(shù)初值的選取對于參數(shù)靈敏度分析而言可以是任意的,但是對于精度適應(yīng)性評價來說則是在參數(shù)率定完成之后,在其優(yōu)值附近進(jìn)行評價;二是評價指標(biāo)的選取在精度適應(yīng)性評價中局限于模擬值與觀測值之間的擬合誤差評價,而靈敏度分析則可以選取任意的模型輸出.

長期以來,地下水?dāng)?shù)值模擬精度評價都是根據(jù)區(qū)域內(nèi)觀測井實測水位與計算水位的差值及擬合程度來判定模型的精度和適應(yīng)性.常用的評價方法是通過最小二乘法計算水位觀測曲線與計算曲線之間的離差平方和來進(jìn)行定量分析,實際也常常通過直觀比較2條曲線的擬合程度來進(jìn)行分析,主觀評定模型擬合的好壞.但是,主觀評定往往受到評定者的主觀喜好和個人經(jīng)驗等因素的制約,不同建模者將會得到不同的模型率定參數(shù),另外通過客觀的最小二乘法或者對數(shù)及加權(quán)求和方法都存在標(biāo)準(zhǔn)化的問題,不同的觀測井,不同時間和區(qū)域,其比較基準(zhǔn)不一樣,因而不能對所得結(jié)果統(tǒng)一進(jìn)行比較.

Nash-Sutcliffe效率系數(shù)E ns及平均相對誤差E mr在徑流過程的相似性和擬合精度評價方面具有常用離差平方和所不具備的標(biāo)準(zhǔn)化特點[3].然而,目前在地下水領(lǐng)域模型模擬精度評價方面尚無采用該類標(biāo)準(zhǔn)的先例.因此,本文以 Ens和 Emr作為評價標(biāo)準(zhǔn),討論地下水?dāng)?shù)值模型模擬結(jié)果的精度適應(yīng)性,并采用Morris方法[4]對模型率定參數(shù)的適應(yīng)性和參數(shù)靈敏度進(jìn)行評價.

1 適應(yīng)性評價方法

1.1 Morris方法

通過靈敏度分析可以評價模型輸出對模型參數(shù)的響應(yīng)程度.靈敏度分析的方法主要有2種:一種為傳統(tǒng)的擾動分析法;另一種為在現(xiàn)代環(huán)境系統(tǒng)中更為有效的區(qū)域靈敏度分析方法[5].這2種方法有不同的應(yīng)用范圍,實施難度也不同[6].擾動分析法中最常用的是Morris方法.Morris方法選取模型中1個變量Pi,其余參數(shù)值固定不變,在修正方法中采用參數(shù)自變量Pi以設(shè)定好的變幅變化,運行模型得到目標(biāo)函數(shù)y(P)=y(P1,P2,P3,…,Pn)的值,用影響值Si判斷參數(shù)變化對輸出值的影響程度[7-9].Si的計算公式為

式中:yn——參數(shù)變化后的輸出值;y0——參數(shù)變化前的輸出值;P0——初始參數(shù);Δi——參數(shù) Pi相對于參數(shù)P0的變幅.靈敏度判別因子S取多個影響值的平均值[8],即

式中n為模型運行次數(shù).

1.2 模擬精度的適應(yīng)性標(biāo)準(zhǔn)

地下水模型中的輸出變量往往是一系列隨時間變化的量,Morris方法輸出變量yi的選取比較困難.這里以Ens和Emr作為輸出標(biāo)準(zhǔn),一方面可以評價數(shù)值模擬精度的適應(yīng)性,另一方面通過Morris方法可以得到模型參數(shù)的靈敏度.E ns和E mr的計算公式為

式中:Hobsi——某時刻觀測水位;Hcali——某時刻計算水位;ˉHobs——觀測水位平均值;ˉHcal——計算水位平均值.

本文定義SN為數(shù)值模擬精度的適應(yīng)性指標(biāo),定義SP為參數(shù)適應(yīng)性判別因子,SP值的大小反映了該參數(shù)在參數(shù)率定結(jié)果附近的綜合靈敏度,S P越大靈敏度越高.S N和S P的計算公式為

式中 Ei為參數(shù)Pi條件下的Ens或 Emr值.

顯然,計算水位與實測水位一致時,效率系數(shù)取得最大值1.一般情況下,效率系數(shù)的變化范圍為0～1.若效率系數(shù)為負(fù)值,則說明模型精度很差.式(3)是目前水文模擬中最常用的目標(biāo)函數(shù)之一.利用該目標(biāo)函數(shù)可以很好地控制模擬過程.但在洪水過程模擬中,可能洪峰預(yù)報精度較高,整個過程將錯位1～2個時段[2].另外,選取Emr協(xié)同評價模型模擬精度.

1.3 數(shù)值模擬精度適應(yīng)性評價

對復(fù)雜非線性系統(tǒng)而言,參數(shù)初值不同,相同變幅引起的系統(tǒng)響應(yīng)也不一樣,并且在同一位置上增大或減小相同幅度的參數(shù)值,對模型的影響也不一定相同[10-13].對建模者來說,最為關(guān)注的是模型參數(shù)在最優(yōu)值附近的靈敏度問題,即模型模擬精度的適應(yīng)性問題.最優(yōu)值附近微小擾動對模型運行結(jié)果的影響是模型模擬精度適應(yīng)性最直接的評價指標(biāo).精度適應(yīng)性判別因子參照文獻(xiàn)[14]中的標(biāo)準(zhǔn)化靈敏度取值,見表1.

表1 精度適應(yīng)性判別因子取值Table 1 Values of sensitivity factors

本研究以離差平方和最小的擬合結(jié)果作為模擬精度適應(yīng)性評價初值,并采用Morris方法,分別以E ns和Emr作為模型輸出來評價模擬精度的適應(yīng)性和參數(shù)的靈敏度,同時以時間序列的觀測水位靈敏度分析結(jié)果作為靈敏度分析的補充依據(jù),討論參數(shù)變化對模型計算水位的影響.

2 計算步驟及結(jié)果

2.1 數(shù)值模型

本文以在貴州普定進(jìn)行的某次抽水試驗數(shù)據(jù)作為研究基礎(chǔ),通過Visual Modflow建立地下水?dāng)?shù)值模型,并采用該數(shù)值模型模擬該次抽水試驗全過程,以抽水井內(nèi)水位動態(tài)資料作為模型驗證和評價資料.該次抽水歷時4h,水位變化主要分為抽水下降段、基本穩(wěn)定段和水位恢復(fù)段3個過程.該次試驗的抽水井為潛水井,最大水位降深為0.93m,相對于含水層厚度來說變化不大,所以選用單層均值非穩(wěn)定流模型進(jìn)行模擬.模型邊界設(shè)置充分大,主要參數(shù)為滲透系數(shù)K和給水度μ.采用Visual Modflow中參數(shù)估計模塊按照觀測值與計算值的離差平方和最小的標(biāo)準(zhǔn)率定模型參數(shù),結(jié)果為K=309 m/d,μ=0.2.觀測水位與計算水位過程如圖1所示,此參數(shù)條件下Ens=0.845,Emr=5.84%.

2.2 模擬精度適應(yīng)性分析

將K=309m/d,μ=0.2作為適應(yīng)性分析的初值,分別對K和μ進(jìn)行增減變化,具體選擇依據(jù)E ns和E mr的變化情況而定,其中參數(shù)K的變幅范圍為-20%～100%,參數(shù) μ的變幅為-50%～100%.Ens,Emr與K,μ變幅的關(guān)系如圖2和圖3所示.

圖1 觀測水位與數(shù)值模型優(yōu)化參數(shù)下的計算水位過程線Fig.1 Comparison of hydrographs of observed and calculated water levels with optimum parameters

圖2 E ns,E mr與 K變幅的關(guān)系Fig.2 Relationship between E ns,E mr and K

圖3 E ns,E mr與 μ變幅的關(guān)系Fig.3 Relationship between E ns,E mr andμ

從圖2可以看出,以離差平方和最小作為判斷依據(jù)的率定結(jié)果與以Ens作為判斷依據(jù)的率定結(jié)果基本一致,而與以平均相對誤差Emr作為判斷依據(jù)的率定結(jié)果分別相差不到5%和10%.這表明K的最優(yōu)值在初值K0～1.05K0之間選取比較合理,μ的最優(yōu)值在μ0～1.1μ0之間選取比較合理.若參數(shù)增大和減小幅度相同,則參數(shù)增大時模型的精度變幅比參數(shù)減小時模型的精度變幅要小.利用式(5)計算得到的SN與K,μ變幅的關(guān)系如圖4所示;利用式(6)計算得到的K和μ的S P值分別為2.38和1.35.

從圖4可以看出:當(dāng)K和μ增大時,SN變化較小;而當(dāng)K和μ減小時,SN變化很大.對比計算得到的SP值與表1,可以發(fā)現(xiàn)K和μ的靈敏度均極高,并且K的靈敏度高于μ.由此可知,參數(shù)優(yōu)化初值確定后,K和μ的微小變化均會給模擬精度帶來很大影響.

圖4 S N與參數(shù)變幅的關(guān)系Fig.4 Variation of adaptability index S N with parameters K andμ

為了進(jìn)行抽水過程的參數(shù)靈敏度分析,利用式(1)計算了影響值Si,其中yi為對應(yīng)于每一時刻的模型計算水位值,y0則是初始參數(shù)對應(yīng)的某一時刻計算水位值.整個時間序列不同變幅K和μ的Si計算結(jié)果如圖5和圖6所示.

圖5 K的Si計算結(jié)果Fig.5 Sensitivity curves of K

圖6 μ的Si計算結(jié)果Fig.6 Sensitivity curves ofμ

比較圖5和圖6可以看出:(a)在模擬期內(nèi),K的靈敏度均高于μ,其值約為后者的2倍.(b)在模擬期的初始階段,井內(nèi)水位變化及其對含水層的影響范圍較小,K與μ的靈敏度均較小;隨著模型計算水位降深的增加,計算水位對參數(shù)靈敏度指標(biāo)值開始增大,初始階段不敏感的參數(shù),模擬后期也具有較高的靈敏度.(c)在模擬后期,水位降深值又開始減小,隨之參數(shù)的靈敏度又有適當(dāng)?shù)慕档?因此,在率定模型參數(shù)時,應(yīng)注意模型高水位降深期間,水位計算值與觀測值的差別.另外,從參數(shù)不同變幅的靈敏度變化曲線來看,K的變幅在-20%～20%之間時,K的靈敏度值較大,但當(dāng)K的變幅大于或等于20%時,K的靈敏度開始降低;μ的變幅在-50%～-5%之間時,μ的靈敏度值較大,而 μ的變幅在-5%～100%之間時,μ的靈敏度開始降低,靈敏度值較小.由此也可以得出在參數(shù)初值取值不同的條件下,其靈敏度分析結(jié)果有所差異的結(jié)論.

3 結(jié) 論

a.基于Morris方法的逐步改變參數(shù)變幅的適應(yīng)性評價方法可以有效地評價優(yōu)化參數(shù)初值附近的模型模擬結(jié)果的精度適應(yīng)性.以Nash-Sutcliffe效率系數(shù)與平均相對誤差作為模擬結(jié)果的精度評價標(biāo)準(zhǔn),所得到的最優(yōu)值存在差異.

b.地下水模型的時間序列靈敏度分析結(jié)果表明,在模擬期初始階段水位降深值較小時,參數(shù)靈敏度較小,并且參數(shù)變幅不同,參數(shù)靈敏度隨時間的變化過程也不同,甚至差異較大.

[1]郝治福,康紹忠.地下水系統(tǒng)數(shù)值模擬的研究現(xiàn)狀和發(fā)展趨勢[J].水利水電科技進(jìn)展,2006,26(1):77-81.(HAO Zhi-fu,KANG Shao-zhong.Current situation and development trend of numerical simulation of groundwater system[J].Advance in Science and Technology of Water Resources,2006,26(1):77-81.(in Chinese))

[2]李慶揚,王能超,易大義.數(shù)值分析[M].4版.北京:清華大學(xué)出版社,2001.

[3]張建云,王國慶.氣候變化對水文水資源影響研究[M].北京:科學(xué)出版社,2007.

[4]CHRISTIAENSK,FEYEN J.Useof sensitivity and uncertainty measures in distributed hydrological modeling with an application to the MIKE SHE model[J].Water Resources Research,2002,38(9):1-15.

[5]徐愛蘭,姚琪,王鵬.基于太湖數(shù)字流域系統(tǒng)的水質(zhì)模型參數(shù)靈敏度分析[J].水利科技與經(jīng)濟(jì),2007,13(1):17-19.(XUAilan,YAOQi,WANG Peng.Sensitivity analysis for parameters of water quality model based on digital valley system of Taihu Basin[J].Water Conservancy Science and Technology and Economy,2007,13(1):17-19.(in Chinese))

[6]王建平,程聲通.軟計算技術(shù)在環(huán)境復(fù)雜模型參數(shù)識別中的應(yīng)用研究[J].系統(tǒng)工程理論與實踐,2006(2):118-126.(WANG Jian-ping,CHENG Shen-tong.Parameter identification of complicated environment model using the soft-computing approach[J].2006(2):118-126.(in Chinese))

[7]FRANCOSA,ELORZA F J.Sensitivity analysis of distributed environmental simulation models:understanding the model behaviour inhydrological studies at the catchment scale[J].Reliability Engineering and System Safty,2003,79(2):205-218.

[8]ZADOR J,ZSELY I G,TURANYI T.Local and global uncertainty analysis of complex chemical kinetic systems[J].Reliability Engineering System Safety,2006,91(10/11):1232-1240.

[9]郝芳華,任希巖,張雪松,等.洛河流域非點源污染負(fù)荷不確定性的影響因素[J].中國環(huán)境科學(xué),2004,24(3):270-274.(HAO Fang-hua,REN Xi-yan,ZHANG Xue-song,et al.Uncertain affecting factor of the non-point source pollution load[J].China Environmental Science,2004,24(3):270-274.(in Chinese))

[10]黃金良,杜鵬飛,何萬謙,等.城市降雨徑流模型的參數(shù)局部靈敏度分析[J].中國環(huán)境科學(xué),2007,27(4):549-553.(HUANG Jin-liang,DU Peng-fei,HE Wan-qian,et al.Local sensitivity analysis for urban rainfall runoff modeling[J].China Environmental Science,2007,27(4):549-553.(in Chinese))

[11]李森,陳家軍,葉慧海,等.地下水流數(shù)值模擬中隨機(jī)因素的靈敏度分析[J].水利學(xué)報,2006,37(8):977-984.(LI Sen,CHEN Jia-jun,YE Hui-hai,et al.Analysis on sensitivity of stochastic factors in numerical simulation of groundwater flow[J].Journal of Hydraulic Engineering,2006,37(8):977-984.(in Chinese))

[12]束龍倉,王茂枚,劉瑞國,等.地下水?dāng)?shù)值模擬中的參數(shù)靈敏度分析[J].河海大學(xué)學(xué)報:自然科學(xué)版,2007,35(5):491-495.(SHULong-cang,WANGMao-mei,LIU Rui-guo,et al.Sensitivity analysis of parameters in numerical simulation of groundwater[J].Journal of Hohai University:Natural Sciences,2007,35(5):491-495.(in Chinese))

[13]SHU Long-cang,LIU Pei-gui,ONGOR B T.Environmental impact assessment using FORM and groundwater system reliability concept:case study Jining,China[J].Environmental Geology,2008,55(3):661-668.

[14]王海龍,余新曉,武思宏,等.SWAT模型靈敏度分析模塊在黃土高原典型流域的應(yīng)用[J].北京林業(yè)大學(xué)學(xué)報,2007,29(增刊2):238-242.(WANG Hai-long,YUXin-xiao,WU Si-hong,et al.Application of sensitivity analysis module of SWATmodel in typical watershed of the Loess Plateau[J].Journal of Beijing Forestry University,2007,29(S2):238-242.(in Chinese))