基于LH-OAT方法的VG模型參數(shù)敏感性分析

秦 萍，王 正，孫兆軍，禹 昭

(1.寧夏大學(xué)環(huán)境工程研究院，銀川 750021；2. 寧夏大學(xué)土木與水利工程學(xué)院，銀川 750021；3. 寧夏大學(xué)新華學(xué)院，銀川 750021；4.寧夏(中阿)旱區(qū)資源評(píng)價(jià)與環(huán)境調(diào)控重點(diǎn)實(shí)驗(yàn)室，銀川 750021)

0 引 言

VG(van Genuchten)模型是進(jìn)行土壤水鹽運(yùn)移研究和土壤水力參數(shù)率定的重要參考依據(jù)，特別是在分析非飽和土壤動(dòng)態(tài)水力特性與討論土壤入滲等過程中，都取得了非常好的擬合效果[1-3]。

但是，VG模型屬于非線性模型，待定參數(shù)較多，且影響VG模型參數(shù)的因素非常繁雜，導(dǎo)致不同地區(qū)不同結(jié)構(gòu)的土壤顯現(xiàn)出的VG模型參數(shù)截然不同[4]。再加上，這些參數(shù)存在的不確定性會(huì)在模型使用過程中產(chǎn)生積累效果，從而影響模型模擬結(jié)果的可信度[5]。因此，必須先率定出可信度較高的模型參數(shù)才能進(jìn)行應(yīng)用，但是在樣本資料較少且不可靠的情況下進(jìn)行率定時(shí)顯得較為困難，對(duì)于這種情況，研究模型參數(shù)的敏感性就具有重要意義。通過模型參數(shù)敏感性的識(shí)別，可以為反演出可信度更高的參數(shù)提供更有力的保障，同時(shí)，也避免了盲目地進(jìn)行參數(shù)的率定。

分析模型參數(shù)的敏感性是為了實(shí)現(xiàn)通過率定最少的參數(shù)同時(shí)得到最佳模擬結(jié)果的目的，通常有局部型和全局型兩類敏感性分析方法[6]。就局部敏感性分析方法來說，其原理是每次選擇一個(gè)或一類參數(shù)作為變量來檢驗(yàn)?zāi)Ｐ蛿M合結(jié)果，而不考慮其他參數(shù)對(duì)其的影響作用，方法容易實(shí)現(xiàn)，但是方法得出的結(jié)果具有一定的片面性和偶然性[7,8]。而全局敏感性分析則是通過改變多個(gè)參數(shù)綜合檢測它們對(duì)模型運(yùn)行結(jié)果的作用程度，用這種方法來分析非線性模型的敏感性具有較大優(yōu)勢。目前，常用的模型參數(shù)全局敏感性分析方法有擾動(dòng)分析法[9]、Monte Carlo法[10]、Latin-Hypercube(LH)模擬法[11]以及LH-OAT(Latin-Hypercube One-factor at-a-time)法[12]等。但在土壤水力參數(shù)模型應(yīng)用中，有關(guān)參數(shù)全局敏感性分析的成果尚不多見。本文針對(duì)土壤水力參數(shù)模型----VG模型，建立基于LH-OAT的全局參數(shù)敏感性分析方法，并利用2017年寧夏中部平原引黃灌區(qū)林地田間數(shù)據(jù)進(jìn)行建模，以開展VG模型中有關(guān)參數(shù)的敏感性分析及模型校正。

1 研究方法

1.1 LH-OAT參數(shù)敏感性分析方法

本文所使用的LH-OAT方法是將由Mckay[13]等提出的LH抽樣法與Morris[14]提出的隨機(jī)OAT敏感性分析方法相結(jié)合產(chǎn)生的。其中，LH抽樣法的基本思想是通過將滿足均勻分布的多維參數(shù)均分為N層，并在保證每個(gè)參數(shù)僅采樣1次的情況下抽樣N次，最后隨機(jī)組合成N個(gè)LH參數(shù)組。但是，通過較少頻次的LH抽樣法雖然能確保敏感性分析的高效性，但是當(dāng)所有參數(shù)擾動(dòng)時(shí)，并不能確定是哪個(gè)參數(shù)的變化引起了輸出值的變化[15]。而OAT敏感性分析方法的基本思想則是通過運(yùn)用控制變量法來對(duì)模型中的每個(gè)參數(shù)進(jìn)行分析，也就是說當(dāng)分析某一參數(shù)時(shí)假定其他參數(shù)保持一定，并對(duì)這一參數(shù)施加一定擾動(dòng)，模型運(yùn)行若干次之后，就可得到每個(gè)參數(shù)的敏感度，但是輸出結(jié)果的擺動(dòng)幅度與其他參數(shù)的取值密切相關(guān)[16]。此外，LH-OAT方法既發(fā)揮了它們各自的優(yōu)點(diǎn)同時(shí)又規(guī)避了它們的缺點(diǎn)，其基本思想是通過對(duì)待分析參數(shù)中的每個(gè)參數(shù)分層并進(jìn)行LH抽樣以得到若干組數(shù)據(jù)，然后對(duì)每組數(shù)據(jù)進(jìn)行敏感性分析，最終得到每組數(shù)據(jù)中每個(gè)參數(shù)的敏感度，最后求取每組數(shù)據(jù)敏感性的平均值。這種方法可以有效地識(shí)別出對(duì)模型結(jié)果產(chǎn)生作用的主要參數(shù)因子，極大地提高了模型參數(shù)的可信度。通常，歸一化后的參數(shù)敏感性用下式計(jì)算[15]：

(1)

式中：S為模型輸出，S=g(G1,G2,…,Gn)；Gi為影響S的參數(shù)(i=1,2,…,n)。

敏感性的各類別如表1所示[17]。

表1 敏感度取值表

Tab.1 Sensitivity values table

分類敏感度指數(shù)敏感性Ⅰ|I|<0.05不敏感Ⅱ0.05≤|I|<0.02一般敏感Ⅲ0.2≤|I|<1.0敏感Ⅳ|I|>1.0極敏感

分析VG模型的m個(gè)參數(shù)的全局敏感性時(shí)，LH-OAT方法的實(shí)現(xiàn)步驟如下。

步驟1 在設(shè)定的m個(gè)參數(shù)值域內(nèi)，采用LH抽樣并隨機(jī)組合生成N個(gè)LH抽樣點(diǎn)；

步驟2 采用OAT方法生成N×(m+1)個(gè)參數(shù)組；

步驟3 針對(duì)步驟3生成的參數(shù)組，重復(fù)運(yùn)行VG計(jì)算主程序N×(m+1)次，輸出變量設(shè)定為土壤含水率θ(h)，土壤相對(duì)飽和度Se，土壤非飽和導(dǎo)水率K(Se)；

步驟4 針對(duì)各變量局部敏感度按照式(2)計(jì)算每個(gè)參數(shù)的相對(duì)敏感性。

Si,j=

j∈[1,n]

(2)

式中：Mj表示為第j個(gè)采樣點(diǎn)對(duì)應(yīng)的目標(biāo)函數(shù)；ei為第i個(gè)參數(shù)；αi為參數(shù)ei的極小變化程度。

步驟5 取參數(shù)ei的最終敏感性指標(biāo)Si為各個(gè)采樣點(diǎn)ei參數(shù)的相對(duì)敏感性平均值，用式(3)計(jì)算。

(3)

LA-OAT方法的計(jì)算流程圖如圖1所示。

i為參數(shù)編號(hào)i=1,2,…,m；j為LH抽樣層編號(hào)j=1,2,…,N圖1 LA-OAT方法流程圖Fig.1 Flow chart of LA-OAT method

1.2 VG模型

描述土壤特征曲線的模型很多，其中VG模型是與實(shí)測值吻合度最高的模型，其模型中的參數(shù)意義明確且受到了廣泛應(yīng)用。Mualem通過對(duì)大量室內(nèi)試驗(yàn)進(jìn)行研究得到土壤飽和導(dǎo)水率Ks與土壤水分特征曲線來估算非飽和導(dǎo)水率模型[18]。van Genuchten將其導(dǎo)出的水分特征曲線形式與Mualem模型相結(jié)合，得出了待定解析形式的VG模型。具體數(shù)學(xué)函數(shù)表達(dá)式為：

(4)

(5)

(6)

式中：θ(h)為土壤體積含水率，cm3/cm3；h為吸力水頭，cm；θr為土壤殘余體積含水率，cm3/cm3；θs為土壤飽和體積含水率，cm3/cm3；α、n是曲線性狀參數(shù)，其中m=1-1/n，α是與土壤物理性質(zhì)有關(guān)的參數(shù)，cm-1；K(Se)為土壤非飽和導(dǎo)水率，cm/min；

Ks為土壤飽和導(dǎo)水率，cm/min；Se為土壤相對(duì)飽和度；l為經(jīng)驗(yàn)擬合參數(shù)，通常取平均值0.5。

1.3 模型中參數(shù)的不確定性

顯然，從式(4)、式(5)、式(6)中可以看出，只要知道參數(shù)θr、θs、α、n、Ks、l就可求出VG模型中的土壤體積含水率、Se為土壤相對(duì)飽和度、K(Se)為土壤非飽和導(dǎo)水率。現(xiàn)在對(duì)上述6個(gè)參數(shù)分別作不確定性的初步分析。研究表明[19]，θr在不同質(zhì)地條件下的數(shù)值較小，而且變化范圍不大，可以利用PTFs法進(jìn)行測定，具有較強(qiáng)的確定性；而θs與θr存在共線關(guān)系，而且θs范圍變化更小，也具有較強(qiáng)的確定性；α、n和Ks具有強(qiáng)烈的非線性關(guān)系，其中，不同土壤物理性質(zhì)的α數(shù)值變化很大，存在很大的不確定性，而n的數(shù)值因土壤質(zhì)地和土壤結(jié)構(gòu)不同而變化，造成了其不確定性，另外Ks的數(shù)值較小，變化范圍較大，而且隨著α的增加，Ks快速增加，具有很強(qiáng)的不確定性。以研究區(qū)4種典型土壤的VG模型參數(shù)為例說明α與Ks之間的關(guān)系，數(shù)據(jù)如表2所示。為了解決α、Ks數(shù)值小且變化范圍大的問題，本文采用了lgα、lgKs及l(fā)gKs/lgα來縮小參數(shù)的范圍； 為一確定的經(jīng)驗(yàn)參數(shù)，通常取0.5。

2 研究區(qū)域概況

本研究區(qū)為寧夏吳忠樹新林場2017年林地田間試驗(yàn)區(qū)(38°36′N，105°56′E)，面積為2.6 htm2。多年平均降水量約為260.7 mm，無霜期176 d，年平均氣溫8.3～8.6 ℃，屬于中溫干旱氣候區(qū)，晝夜溫差大，四季分明。試驗(yàn)地土壤類型為粗質(zhì)灰鈣土，土壤質(zhì)地以砂壤土為主，有少數(shù)砂土和壤土，土壤密度為1.48～1.67 g/cm3，有機(jī)質(zhì)含量較高，而深層土壤則以黏性土壤為主。試驗(yàn)時(shí)間為2017年11月7-27日，開展了土壤含水率和土壤物理性質(zhì)的測定。試驗(yàn)內(nèi)容主要包括：①土壤物理性質(zhì)：根據(jù)試驗(yàn)地土壤采樣點(diǎn)的分層取樣樣本，利用激光粒度儀分析試驗(yàn)地0～100 cm深度的土壤機(jī)械組成；②土壤含水率：采用烘干法分5層測定(0～20、20～40、40～60、60～80、80～100 cm)，監(jiān)測周期為6～8 d。

表2 研究區(qū)4種典型土壤的VG模型參數(shù)

Tab.2 Soil hydraulic parameters of the VG model for 4 soils

土壤類型θr/(cm3·cm-3)θs/ (cm3·cm-3)α/cm-1nKs/(cm·min-1) lgαlgKslgKs/lga砂土0.045 00.432 10.145 22.682 10.495 0-0.838 6-0.305 42.745 9砂質(zhì)黏壤土0.100 00.398 40.124 52.281 40.243 2-0.906 5-0.614 41.475 4砂壤土0.065 00.435 20.059 41.487 50.021 8-1.226 2-1.661 50.738 0壤土0.078 00.467 30.036 31.566 50.017 3-1.443 7-1.761 20.819 7

模型設(shè)置時(shí)，土壤剖面厚度設(shè)定為100 cm，依據(jù)剖面土壤質(zhì)地分為0～60 cm(砂土層)和60～100 cm(黏土層)的2個(gè)水平分層，記為a、b層，初始條件為模擬初始日期的實(shí)測土壤含水率。

3 分析與討論

VG模型需要率定的參數(shù)有θr、θs、α、n、Ks、l等。本研究VG模型取值依照RETC軟件進(jìn)行確定，各個(gè)參數(shù)的取值范圍如表3所示。

表3 VG模型參數(shù)取值范圍

Tab.3 The range of the parameters of the VG model

土壤類型θr/(cm3·cm-3)θs/(cm3·cm-3)α/cm-1nKs/(cm·min-1)砂土0.045 0～0.057 00.428 5～0.430 30.036 3～0.145 21.890 2～2.683 20.073 6～0.495 0砂質(zhì)黏壤土0.068 0～0.100 00.378 7～0.381 20.008 0～0.390 01.091 5～1.482 50.003 3～0.021 4

本文針對(duì)VG模型不同的輸出結(jié)果(θ(h)、Se、K(S3))，分別計(jì)算土壤a層、b層的θr、θs、α、n、Ks、l等參數(shù)，共獲得了11個(gè)參數(shù)的敏感性排序結(jié)果，如表4所示。

按照表1給出的敏感性取值，極敏感參數(shù)，排序?yàn)?；敏感參數(shù)，排序?yàn)?～6；一般敏感參數(shù)，排序?yàn)?～11；不敏感參數(shù)，排序?yàn)?2。取出各參數(shù)排序的最小值作為全局敏感性分析的排序結(jié)果，如表5所示。

表4 參數(shù)敏感性順序

Tab.4 The order of parameters sensitivity

序號(hào)123456 參數(shù)θsa、θsbθra、θrbαa、αbKsa、Ksbna、nbl

注：下標(biāo)a、b表示土層

表5 參數(shù)敏感性分析結(jié)果

Tab.5 The results of sensitivity analysis

參數(shù)敏感性等級(jí)S3K(S3)θd1(0～20 cm)θd2(20～40 cm)θd3(40～60 cm)θd4(60～80 cm)θd5(80～100 cm)θsa11111211αa13224122na22443333Ksa2113354116αb251179545θra381066964Ksb510810111187nb6771187910θsb86697879θrb99981010108l1212121212121212

由于模型中各個(gè)參數(shù)參與的過程不同，對(duì)于不同輸出變量的敏感性也有較大差異，全局敏感度的柱狀圖如圖2所示。

圖2 全局敏感度 Fig.2 The values of global sensitivity

另外，分別計(jì)算各個(gè)參數(shù)對(duì)不同輸出變量敏感度的平均值作為該參數(shù)的全局敏感度，如圖3所示，相比較而言，極敏感參數(shù)的敏感度指數(shù)明顯大于敏感參數(shù)和一般敏感參數(shù)，與表1所指出的數(shù)據(jù)范圍相一致，而不敏感參數(shù)的全局敏感度為0。以Se為輸出變量時(shí)，其對(duì)θsa、αa和na呈高度的敏感性，而對(duì)αb、θra、nb和θsb表現(xiàn)出較高的敏感性，對(duì)Ksa、Ksb和l表現(xiàn)為一般敏感或不敏感。以K(Se)為輸出變量時(shí)，其對(duì)θsa表現(xiàn)最為敏感，對(duì)αa、na、Ksa表現(xiàn)出較高的敏感性，對(duì)其他參數(shù)敏感性較弱。

圖3 輸出為Se的全局敏感度Fig.3 The values of global sensitivity

圖4 輸出為K(Se)的全局敏感度Fig.4 The values of global sensitivity

下面分析輸出為土壤含水量時(shí)的模型全局敏感度，結(jié)果如圖5所示。本研究的試驗(yàn)區(qū)采用漫灌方式，水分充足導(dǎo)致土壤分層含水量主要與土壤水力參數(shù)相關(guān)。考慮到受土壤蒸發(fā)和植物需水等因素影響，再加上土壤樣本的監(jiān)測點(diǎn)主要集中在a層，所以土壤含水量對(duì)a層的土壤水力參數(shù)均呈現(xiàn)出較強(qiáng)的敏感性，土壤水分特征曲線的主要參數(shù)θsa、αa和na表現(xiàn)出強(qiáng)烈的敏感性，同時(shí)b層土壤部分參數(shù)也對(duì)a層土壤含水量輸出結(jié)果呈現(xiàn)出了一定的敏感性，這是因?yàn)榈叵滤挠潟?huì)導(dǎo)致底層土壤產(chǎn)生滲漏和水分補(bǔ)充等過程，從而導(dǎo)致了底層土壤水力參數(shù)對(duì)上層土壤的含水量產(chǎn)生作用。

圖5 不同深度土層含水量的全局敏感性Fig.5 Global sensitivity of soil water content at different depths

4 模型校驗(yàn)

為了得到更好地?cái)M合效果，對(duì)上述敏感參數(shù)實(shí)施手動(dòng)修正，以期達(dá)到最佳擬合效果。經(jīng)過多次調(diào)整，修正后的敏感參數(shù)值如表6所示。將修正后的敏感參數(shù)與極敏感參數(shù)一起應(yīng)用于VG模型中，模型模擬效果采用平均相對(duì)誤差(MRE)、均方差(RMSE)及決定系數(shù)(R2)來評(píng)價(jià)[19]。某樣本點(diǎn)修正后的土壤含水率實(shí)測值與模擬值對(duì)比如圖6所示。

通過對(duì)比后容易得知，不同監(jiān)測時(shí)期不同土層深度的土壤含水率有顯著變化，且模擬值較好地反映了實(shí)測值的變化趨勢，吻合度較高。由于第1次采集土樣發(fā)生在冬灌后的第7 d后，且上層土壤受蒸發(fā)、下滲和植物吸水等因素的影響，所以土壤下層的含水率要高于上層。隨時(shí)間的推移，各層的土壤含水率均有所下降，這也與實(shí)際相符合。

表6 VG模型參數(shù)對(duì)比

Tab.6 The range of the parameters of the VG model

土壤類型θra/(cm3·cm-3)θsa/ (cm3·cm-3)θsb/ (cm3·cm-3)αa/cm-1αb/cm-1nanb修正前0.0450.436 80.379 20.040 80.038 61.3791.400修正后0.0510.509 80.382 40.042 60.047 31.4061.366

圖6 土壤樣本含水率實(shí)測值與模擬值對(duì)比 Fig.6 Comparison between the simulated and the measured value

表7 VG模型檢驗(yàn)指標(biāo)

Tab.7 The test indexes of VG model

土壤類型MRE/%RMSER2Se2.648 70.196 80.924 3K(Se)3.023 50.236 50.864 4土壤含水率1.960 70.046 60.974 4

由表7的VG模型檢驗(yàn)指標(biāo)可以看出，MRE、RMSE都較小，而且R2接近于1，較為理想，因此模擬結(jié)果可以接受。結(jié)合土壤樣本含水率實(shí)測值與模擬值對(duì)比(圖6)、模擬檢驗(yàn)指標(biāo)(表7)可知，修正后的敏感參數(shù)應(yīng)用于VG模型能較好地模擬該研究區(qū)內(nèi)的不同土層含水率。

5 結(jié) 語

本研究將LH-OAT方法與VG模型相結(jié)合，討論了適用于VG模型參數(shù)全局敏感性分析方法，并將其應(yīng)用于寧夏引黃灌區(qū)2017年的土壤含水率數(shù)據(jù)進(jìn)行了方法檢測和參數(shù)敏感性分析，并對(duì)分析后的敏感參數(shù)進(jìn)行了修正。結(jié)論如下：

(1)VG模型中參數(shù)θsa、αa和na表現(xiàn)極為敏感，而αb、θra、nb和θsb較為敏感，Ksa、Ksb和l表現(xiàn)為一般敏感或不敏感。

(2)定量討論了參數(shù)對(duì)各輸出變量(θ(h)、Se、K(Se))的敏感性程度，對(duì)不同輸出呈現(xiàn)顯著的差異。

(3)根據(jù)全局敏感性分析的結(jié)果，修正敏感參數(shù)后對(duì)模型進(jìn)行了校驗(yàn)，校驗(yàn)的結(jié)果較為理想。該方法具有很強(qiáng)的普適性，為提高后續(xù)模型參數(shù)率定效率提供了依據(jù)。