基于靈敏度分析的地下水?dāng)?shù)值模擬精度適應(yīng)性評(píng)價(jià)

魯程鵬,束龍倉(cāng),劉麗紅,劉佩貴,董貴明

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

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

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

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

1 適應(yīng)性評(píng)價(jià)方法

1.1 Morris方法

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

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

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

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

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

式中:Hobsi——某時(shí)刻觀測(cè)水位;Hcali——某時(shí)刻計(jì)算水位;ˉHobs——觀測(cè)水位平均值;ˉHcal——計(jì)算水位平均值.

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

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

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

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

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

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

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

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

2.1 數(shù)值模型

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

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

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

圖1 觀測(cè)水位與數(shù)值模型優(yōu)化參數(shù)下的計(jì)算水位過程線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é)果基本一致,而與以平均相對(duì)誤差Emr作為判斷依據(jù)的率定結(jié)果分別相差不到5%和10%.這表明K的最優(yōu)值在初值K0～1.05K0之間選取比較合理,μ的最優(yōu)值在μ0～1.1μ0之間選取比較合理.若參數(shù)增大和減小幅度相同,則參數(shù)增大時(shí)模型的精度變幅比參數(shù)減小時(shí)模型的精度變幅要小.利用式(5)計(jì)算得到的SN與K,μ變幅的關(guān)系如圖4所示;利用式(6)計(jì)算得到的K和μ的S P值分別為2.38和1.35.

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

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

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

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

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

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

3 結(jié) 論

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

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

[1]郝治福,康紹忠.地下水系統(tǒng)數(shù)值模擬的研究現(xiàn)狀和發(fā)展趨勢(shì)[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]李慶揚(yáng),王能超,易大義.數(shù)值分析[M].4版.北京:清華大學(xué)出版社,2001.

[3]張建云,王國(guó)慶.氣候變化對(duì)水文水資源影響研究[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]王建平,程聲通.軟計(jì)算技術(shù)在環(huán)境復(fù)雜模型參數(shù)識(shí)別中的應(yīng)用研究[J].系統(tǒng)工程理論與實(shí)踐,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]郝芳華,任希巖,張雪松,等.洛河流域非點(diǎn)源污染負(fù)荷不確定性的影響因素[J].中國(guó)環(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].中國(guó)環(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é)報(bào),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]束龍倉(cāng),王茂枚,劉瑞國(guó),等.地下水?dāng)?shù)值模擬中的參數(shù)靈敏度分析[J].河海大學(xué)學(xué)報(bào):自然科學(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é)報(bào),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))