桂晗亮, 張春萍,, 武治國, 張家銓

(1.武漢工程大學(xué) 光電信息與能源工程學(xué)院, 湖北 武漢 430000; 2.武漢新烽光電股份有限公司, 湖北 武漢 430000)

近些年來，隨著城市工業(yè)化進(jìn)程的加快，環(huán)境受到人類活動(dòng)影響的程度不斷加劇，其中水環(huán)境問題日漸嚴(yán)重，人類在享受工業(yè)化所帶來的效益同時(shí)，不得不開始重視水環(huán)境問題。流域水文模型具有成本低、便于情景分析等優(yōu)點(diǎn)，現(xiàn)被廣泛用于水文分析計(jì)算、生態(tài)保護(hù)等方面[1]。而流域水文模型可按對(duì)流域水文過程描述的離散程度分為：集總式、分布式和半分布式[2-3]。由于流域模擬問題的復(fù)雜性等原因，使得半分布式模型在實(shí)踐的過程中具有更廣闊的應(yīng)用前景和更大的推廣價(jià)值。而HSPF(hydrological simulation program-Fortran)模型作為半分布式模型的經(jīng)典代表，國外對(duì)于HSPF模型的應(yīng)用已經(jīng)非常廣泛，研究者[4,9,13,23-24]利用HSPF模型對(duì)降雨徑流為主的水文過程、氣候變化以及土地利用變化對(duì)水文過程的影響開展了很多研究。Ourania Tzoraki等[5]利用HSPF模型模擬喀斯特流域的水文過程。Diaz-Ramirez等[6]的研究則表明，HSPF能夠適用于熱帶島嶼河流流域水文模擬。Luo C.等[24]通過HSPF模型模擬土地利用的變化對(duì)Xitiaoxi流域水文水質(zhì)的影響。然而HSPF模型在國內(nèi)的關(guān)注和應(yīng)用卻很少，近年來才陸續(xù)有研究將該模型應(yīng)用于丘陵地區(qū)太湖流域[7]、半干旱地區(qū)媯水河流域等[8]，但鮮有研究證實(shí)該模型在熱帶沿海流域的適用性。此外，也有學(xué)者對(duì)于HSPF模型的不確定性問題進(jìn)行了深入研究，Igor Iskra等[9]利用蒙特卡洛—拉丁超立方(MC-LHS)、響應(yīng)面和矩估計(jì)3種方法對(duì)HSPF進(jìn)行了不確定性分析。程曉光等[8]利用GLUE-MC方法證實(shí)了HSPF模型不確定性與徑流量有關(guān)。而龐樹江等[10]利用CLUE分析HSPF模型的不確定性，并提出不確定性大小存在季節(jié)差異性。但是，對(duì)于降雨量影響模型不確定性的研究卻很少見。

因此，為進(jìn)一步探究HSPF模型在不同地區(qū)的適用性和不確定性的影響因素。本研究通過建立熱帶沿海地區(qū)三亞河流域的HSPF水文模型，選取2017—2019年3 a的徑流量進(jìn)行模擬，探究HSPF模型的敏感度以及不同降雨量條件下模型的不確定性，以期為HSPF模型在不同流域和地區(qū)的應(yīng)用提供參考和依據(jù)。

1 研究區(qū)域與方法

1.1 研究區(qū)域概述

三亞河位于三亞市南部，由六羅水、水蛟溪、半嶺水三條河組成，以六羅水為主流。三亞河由北向南穿城而過，經(jīng)由三亞港入海。地理位置位為109°28′—109°31′E，18°14′—18°22′N，總流程28.8 km，流域面積337 km2。氣候?yàn)闊釒ШＱ笮约撅L(fēng)氣候，年平均溫度為25.7 ℃，年平均降雨量為1 347.5 mm，且降雨主要集中在7—11月。土地利用類型包括林地、農(nóng)業(yè)用地、城市用地和草地。

1.2 模型架構(gòu)

建立水文模型輸入的數(shù)據(jù)包括時(shí)間序列WDM數(shù)據(jù)(氣象、水文數(shù)據(jù))和GIS數(shù)據(jù)(地形DEM、土地利用數(shù)據(jù))(詳見表1)。利用BASINS模型所嵌套的WDMUtil軟件將氣象和水文數(shù)據(jù)整合成時(shí)間序列WDM文件，同時(shí)通過BASINS和GIS軟件分別進(jìn)行流域的劃分(泰森多邊形法)和邊界的提取，將三亞河流域劃分為28個(gè)子流域，HSPF建模需要的地理數(shù)據(jù)，由BASINS向WinHSPF轉(zhuǎn)化，進(jìn)而生成WinHSPF的UCI運(yùn)行文件，并利用PEST軟件對(duì)模型參數(shù)進(jìn)行尋優(yōu)。

表1 三亞河流域建模數(shù)據(jù)內(nèi)容

模型的率定和驗(yàn)證是水文預(yù)報(bào)模型中必不可少的部分。研究采用納什系數(shù)(NSE)和相對(duì)誤差(Re)對(duì)模型擬合程度進(jìn)行評(píng)價(jià)，NSE越接近1，Re越接近0，說明模型擬合的程度就越好。

(1)

(2)

1.3 參數(shù)敏感度分析

參數(shù)敏感度分析是研究參數(shù)變化所引起的模型響應(yīng)，是模型不確定性分析的重要內(nèi)容[11]，也是進(jìn)行不確定性分析的前提工作之一。而HSPF模型水文模塊的參數(shù)眾多，且大部分參數(shù)具有物理背景，因此為提高不確定性研究的效率，需要先通過參數(shù)的敏感度分析，篩選出水文模塊中對(duì)模擬有影響的參數(shù)。基于擾動(dòng)分析法的Morris篩選法是目前應(yīng)用較廣的一種敏感度分析方法[12]。本研究采用修正的Morris篩選法，自變量以固定步長變化(±20%，±15%，±10%，±5%)，運(yùn)用公式(3)計(jì)算出Morris系數(shù)的平均值用以評(píng)價(jià)參數(shù)的敏感度大小。

(3)

式中：S為參數(shù)相對(duì)敏感度(0≤|S|<0.05,不敏感； 0.05≤|S|<0.2,一般敏感； 0.2≤|S|<1,敏感； |S|≥1,極敏感)；Yi，Yi+1為模型第i,i+1次運(yùn)行輸出值；Y0為參數(shù)率定后計(jì)算結(jié)果初始值；Pi，Pi+1為第i，i+1次模型運(yùn)算參數(shù)值相對(duì)于率定后參數(shù)值的變化百分率；n為模型運(yùn)行次數(shù)。

HSPF模型包括3個(gè)主要模塊：PERLND(透水區(qū))，IMPLND(不透水區(qū))和RCHRES(河道)，其中IMPLND模塊僅分析地表徑流，而PERLND模塊主要分析地表徑流、壤中流和地下水水文過程。模型自上而下分為樹冠層、植被層、各土壤層(包括表層土壤、上土壤層、下土壤層、地下水涵養(yǎng)層)，降雨在6個(gè)垂直的存儲(chǔ)層進(jìn)行分配，最終進(jìn)入河道[7]。研究主要通過HSPF的水文模塊模擬三亞河流域的徑流量，結(jié)合相關(guān)研究[11,13-15]，選取對(duì)HSPF水文模塊有影響的8個(gè)參數(shù)及其取值范圍(表2)，利用修正的Morris篩選法，計(jì)算各參數(shù)在該流域的相對(duì)敏感度。

表2 對(duì)HSPF水文模塊有影響的參數(shù)敏感度分析

1.4 MC-LHS分析方法

MC-LHS方法是將蒙特卡洛(Monte-Carlo, MC)模擬和拉丁超立方抽樣(Latin hypercube sampling, LHS)結(jié)合的一種不確定性分析方法。蒙特卡洛方法則是一種以概率統(tǒng)計(jì)理論為基礎(chǔ)的數(shù)值分析方法，常用于解決確定性和不確定性問題[3]。蒙特卡洛模擬的主要步驟分為： ①對(duì)輸入?yún)?shù)范圍內(nèi)隨機(jī)抽樣； ②將隨機(jī)抽樣結(jié)果代入模型運(yùn)行； ③對(duì)運(yùn)行結(jié)果進(jìn)行統(tǒng)計(jì)分析。但是，蒙特卡洛模擬中的隨機(jī)抽樣往往不能保證參數(shù)的抽樣點(diǎn)服從均勻分布，且需要非常龐大的樣本量，而基于分層隨機(jī)抽樣的拉丁超立方抽樣能夠使樣本的覆蓋率更好，同時(shí)能減少抽樣次數(shù)[16]。因此將拉丁超立方抽樣與蒙特卡洛結(jié)合起來的MC-LHS方法用于模型不確定性分析的效率更高。本研究以納什系數(shù)為目標(biāo)函數(shù)，選取敏感度分析結(jié)果中的參數(shù)，采用拉丁超立方抽樣取得500個(gè)樣本，代入HSPF模型中運(yùn)行并計(jì)算納什系數(shù)，通過公式(4)—(5)對(duì)目標(biāo)值納什系數(shù)進(jìn)行區(qū)間估計(jì)和變異系數(shù)計(jì)算。

(4)

(5)

式中：X為目標(biāo)值；E(X)，D(X)分別為隨機(jī)變量X的數(shù)學(xué)期望和方差；ε為任意給定的正數(shù)；Cv為變異系數(shù)； SD， MN分別表示隨機(jī)變量X的標(biāo)準(zhǔn)偏差和均值。

2 結(jié)果與分析

2.1 模型率定與驗(yàn)證

研究選用三亞河流域2017—2019年日徑流量為模擬對(duì)象，其中2017年1月1日至2018年12月31日為率定期，2019年1月1日至2019年12月31日為驗(yàn)證期，選用納什系數(shù)和相對(duì)誤差來評(píng)價(jià)模型擬合程度，結(jié)果詳見表3。結(jié)合圖1的模型擬合過程線，率定期和驗(yàn)證期的納什系數(shù)分別為0.93和0.98，相對(duì)誤差分別為0.87%和0.21%，就水文模型模擬而言，一般認(rèn)為納什系數(shù)大于0.90和相對(duì)誤差小于10%就表明模型模擬的效果很好[17]。從率定期和驗(yàn)證期模型擬合的結(jié)果來看，建立的HSPF模型對(duì)于日時(shí)間步長的徑流量擬合程度很好。為進(jìn)一步驗(yàn)證模型的適用性，將模型用于月時(shí)間步長的徑流量的模擬，如表3和圖1所示。結(jié)果說明HSPF模型適用于不同時(shí)間步長沿海熱帶流域的模擬，能夠較為真實(shí)地反映實(shí)際水文過程。

表3 HSPF模型率定和驗(yàn)證的模擬結(jié)果

圖1 不同時(shí)間步長降雨徑流過程線及散點(diǎn)圖

2.2 參數(shù)敏感度分析結(jié)果

以模型擬合程度指標(biāo)NSE為目標(biāo)值，利用修正的Morris篩選法計(jì)算參數(shù)敏感度(表2)。分析表2可知，參數(shù)敏感度依次為：AGWRC，LZSN，INFILT，INTWF，UZSN，IRC，LZETP和DEEPFR，其中相對(duì)敏感度較大的參數(shù)為AGWRC，LZSN，INFILT和INTWF，AGWRC的相對(duì)敏感度為-0.69(敏感)，LZSN的相對(duì)敏感度為-0.05(一般敏感)，INFILT和INTWF的相對(duì)敏感度為-0.04，余下參數(shù)的相對(duì)敏感度較低。

AGWRC是地下水消退系數(shù)，控制地下水的退水過程進(jìn)而影響產(chǎn)流量；LZSN是下土壤層的額定蓄積，主要影響河道入流量[18]。HSPF降雨匯流過程主要是由地表徑流、壤中流和基流組成，熱帶流域獨(dú)特的氣候條件有可能影響土壤含水量、下滲率等特性，使得地表徑流和壤中流在降雨匯流過程中的作用減弱，此時(shí)水文過程主要受基流影響。

HSPF模型的水文循環(huán)包括植被截留、地表徑流、土壤水分配、地下水運(yùn)動(dòng)、蒸散發(fā)等部分。通過與國內(nèi)外對(duì)HSPF模型水文模塊參數(shù)敏感度分析的研究結(jié)論進(jìn)行對(duì)比(表4)，發(fā)現(xiàn)HSPF模型參數(shù)敏感度的相對(duì)大小具有空間差異性，不同流域的敏感參數(shù)和相對(duì)敏感度大小并不完全相同，本研究中較為敏感的參數(shù)依次為AGWRC，LZSN，INFILT和INTWF，而其他研究[7,11]中敏感的參數(shù)更多。結(jié)合模型的水文機(jī)理分析，模型在不同流域的模擬中存在差異性[19]，可能與各流域不同的地形地貌、土地利用等條件有很大關(guān)系。例如，將本研究與半干旱半濕潤的媯水河流域[8]進(jìn)行對(duì)比，發(fā)現(xiàn)在參數(shù)的相對(duì)敏感度大小上就存在著很大的差異性，媯水河流域地處北溫帶，屬于大陸性季風(fēng)氣候，地形多以丘陵為主，本研究中的熱帶沿海流域?qū)儆跓釒ШＱ蠹撅L(fēng)性氣候，多以城鎮(zhèn)、山地為主，在地形地貌、氣候、城鎮(zhèn)化比例等方面明顯的差異，能夠影響地表徑流和土壤層中的再分配過程。此外也有研究指出流域特征和氣候特征等相似的流域之間，敏感性參數(shù)類似[20]。

表4 國內(nèi)外HSPF參數(shù)敏感度分析對(duì)比

2.3 不同降雨量對(duì)模型不確定性的影響

選取敏感參數(shù)進(jìn)行拉丁超立方抽樣并代入HSPF模型中運(yùn)行，依據(jù)日降雨量劃分成小雨(0～10 mm)、中雨(10～25 mm)、大雨(25～50 mm)和暴雨(50 mm以上)4種降雨類型，以模型擬合程度評(píng)價(jià)指標(biāo)納什系數(shù)為目標(biāo)函數(shù)用于不確定分析，分別計(jì)算離散系數(shù)和置信水平為80%，90%下模型區(qū)間(表5)。由表5可知，4種類型降雨的Cv值分別為大雨(-11.268)、暴雨(-4.074)、中雨(0.104)、小雨(0.102)，說明HSPF模型對(duì)于不同降雨量模擬的不確定性存在差異，降雨量小的時(shí)段(≤25 mm)模型的擬合程度相比較降雨量大的時(shí)段(>25 mm)模型的擬合程度更好，且降雨量越大，不確定性越大，模型就越不穩(wěn)定。分別對(duì)比置信水平為80%和90%下目標(biāo)值的變化區(qū)間，模型區(qū)間也與降雨量有關(guān)，降雨量較小時(shí)(≤25 mm)不同置信水平下的變化區(qū)間遠(yuǎn)小于降雨量大時(shí)(>25 mm)。

表5 不同降雨類型HSPF不確定性計(jì)算結(jié)果

降雨量不同時(shí)，模型不確定性大小和區(qū)間有很大差異(圖2)。圖2表明HSPF模型的不確定性大小和降雨量有關(guān)。之所以出現(xiàn)這種情況，與HSPF模型的產(chǎn)流機(jī)制有關(guān)，模型中將透水地段水量平衡過程概化成凈雨到達(dá)地表形成地表徑流，產(chǎn)生入滲并實(shí)現(xiàn)各土壤層中的再分配，不透水地段水文過程為凈雨產(chǎn)生地表徑流的過程；降雨經(jīng)地表植被和洼地截留后產(chǎn)生地表徑流、壤中流和基流，落地雨在土壤中的再分配過程與雨強(qiáng)等因素密切相關(guān)。當(dāng)雨量較小時(shí)，降雨大部分直接經(jīng)過地表植被等截留后直接形成了地表徑流，而少部分或沒有凈雨進(jìn)入土壤中再分配，這時(shí)模型產(chǎn)流主要是地表徑流，但是雨量較大時(shí)，模型的產(chǎn)流包括地表徑流、壤中流和基流。因此，降雨量可用于解釋其他研究[10]提出的HSPF模型不確定性的季節(jié)差異，同樣在其他水文預(yù)報(bào)模型的研究中也證實(shí)了降雨量是影響水文模型不確定性的因素之一[21-22]。

3 結(jié) 論

(1) 通過HSPF模型對(duì)三亞河流域的日流量進(jìn)行模擬，分別以2017和2018年的日徑流量為率定期，2019年為驗(yàn)證期，率定期和驗(yàn)證期的納什系數(shù)為0.93和0.98，相對(duì)誤差為0.87%和0.21%。結(jié)果表明所建立的HSPF模型能夠很好的模擬三亞河流域?qū)嶋H的水文過程。此外，進(jìn)一步驗(yàn)證了HSPF模型對(duì)月時(shí)間步長的模擬，證實(shí)了在沿海熱帶流域的適用性。

圖2 不同程度降雨目標(biāo)值箱線圖

(2) 在HSPF模型完成率定的基礎(chǔ)上，利用Morris篩選法對(duì)模型水文模塊參數(shù)進(jìn)行敏感度分析，三亞河流域?qū)搅鬟^程影響的敏感參數(shù)依次為：AGWRC，LZSN，INFILT，INTWF，UZSN，IRC，LZETP，DEEPFR，其中AGWRC最敏感，而LZETP和DEEPFR對(duì)該流域徑流模擬幾乎無影響；同時(shí)與國內(nèi)外研究中的結(jié)論進(jìn)行對(duì)比，發(fā)現(xiàn)參數(shù)敏感度在不同流域不同。

(3) 利用MC-LHS方法進(jìn)行不確定分析，發(fā)現(xiàn)模型模擬的不確定性與降雨量之間相關(guān)性明顯。降雨量越大，模型模擬的置信區(qū)間和變異系數(shù)Cv越大，模型的不確定性就越大，說明模型就越不穩(wěn)定。同時(shí)證實(shí)了降雨量是影響HSPF模型不確定性的因素之一，進(jìn)一步解釋降雨量的變化是造成HSPF模型不確定存在季節(jié)差異性的原因之一。