陽 揚，翟祿新，賈艷紅，薛開元

（1.廣西師范大學環(huán)境與資源學院，廣西桂林 541004；2.珍稀瀕危動植物生態(tài)與環(huán)境保護教育部重點實驗室，廣西桂林 541004）

0 引 言

氣候和下墊面變化是影響流域水文循環(huán)過程改變的主要因素，其中徑流變化對不同時空尺度下氣候和下墊面變化的響應研究是水資源規(guī)劃及其管理的基礎［1，2］。對于描述自然流域的水文模型，模型參數大多反映了流域土壤和植被等的水文物理特性［3-5］，與流域的土壤、植被、地形地質條件等相關。由于以往假定流域特征是處于穩(wěn)定狀態(tài)，即土壤、地形、地貌等下墊面在長時間尺度上是不變的，因而通常率定后的模型參數隨時間固定不變。但實際情況是，人類活動通過水利工程建設、改變土地利用和植被覆蓋等方式直接或間接影響流域對大氣降水的分配以及蒸散發(fā)等，從而影響地表徑流和地下徑流的形成過程，表現(xiàn)出非一致性特征［6-8］。因此，水文學中傳統(tǒng)的靜態(tài)參數假設已不能反映變化環(huán)境下流域的特征，水文模型參數應隨流域特征產生動態(tài)變化。Wagener T 在探討如何模擬環(huán)境變化對水文過程的影響時，提出模型參數與流域特征具有復雜的相關關系，應當是流域多種特征的函數，隨著流域特征的變化而變化［9］。Merz R 等對奧地利273 個流域的水文狀況進行模擬，發(fā)現(xiàn)降雪融雪過程的相關參數、土壤特性參數和產流過程相關參數等均隨時間發(fā)生顯著變化［10］。由此看來，土地利用和植被覆蓋與模型參數的動態(tài)變化存在相關性，即提供一種思路：基于土地利用指標與模型參數的相關性，可以分析人類活動影響下流域特征變化對降水量-徑流機制的影響，從而實現(xiàn)流域變化的動態(tài)定量感知。孫瑜對漢江流域和潮河流域開展研究，認為不同模型的同類參數與環(huán)境要素呈現(xiàn)出同向的相關關系，如蒸發(fā)過程參數均與植被覆蓋條件正相關，產流過程參數均與植被覆蓋條件負相關［4］。熊立華［11］等考慮植被和人類活動影響的水文模型參數時變特征，認為最大蓄水能力Sc與流域植被覆蓋條件、土地利用變化、水利工程修筑等具有較強的相關關系。

本文選用集總式水文模型MODHYDROLOG 對漓江流域徑流開展模擬。該模型結構簡單清晰、便于率定和計算，而且由于集總式水文模型對基礎資料要求較低，可在研究區(qū)基礎地理數據并非特別豐富的情況下應用，降低了開展水文模擬的難度［12-18］。同時該模型考慮了超滲產流和蓄滿產流兩種產流機制，因此在干旱、濕潤流域都能得以廣泛應用。據此，本文在不同時段人類活動和氣候變化情景下，通過時變參數分析流域降水-徑流特征的變化情況，探索劇烈人類活動影響下流域下墊面變化動態(tài)定量感知，為流域非穩(wěn)態(tài)背景下的水文模擬提供新的思路，以提高定量分析環(huán)境變化對于流域水文變量影響的精度，揭示變化環(huán)境下整個流域氣候變化、人類活動、土壤和植被等變化對河川徑流的影響，從而為流域防洪減災、流域生態(tài)系統(tǒng)的穩(wěn)定性評價、水源涵養(yǎng)功能提升等提供科學依據。

1 研究區(qū)概況

漓江起源于貓兒山，流經興安縣、靈川縣、臨桂縣、桂林市區(qū)、陽朔縣、平樂縣城之后與恭城河交匯，流域全長約214 km，面積約為1.23 萬km2，地理坐標為24°38′10″ N～25°53′59″ N、110°07′39″E～110°42′57″E（圖1）［19］。整體地勢北高南低，以漓江為中軸線，呈西北東南向分布，北部和東部為碎屑巖中低山地貌，中南部主要為低山、丘陵、巖溶地貌，具有典型的喀斯特地貌特征［20，21］。本文研究區(qū)為桂林水文站以上流域范圍，集水面積約為2 762 km2。漓江桂林水文站以上流域水利工程較多，已建成的水庫主要有：青獅潭水庫、川江水庫、小溶江水庫、斧子口水庫等［19］。

圖1 漓江流域范圍Fig.1 Scope of Lijiang River

2 數據與方法

2.1 數據收集

本文主要收集桂林水文局的桂林水文站2007-2016年日尺度徑流數據及基本氣象數據。研究還收集了流域內的圖形數據、DEM 數據以及相關的規(guī)劃圖等，主要有2020年桂林市行政區(qū)劃圖、漓江流域30 m 分辨率數字高程，使用2008、2010、2013以及2015年共4年桂林市1∶1 萬土地利用現(xiàn)狀圖，其他年份采用年內線性插值方法推求獲取。土地利用分類參照表1。

表1 中國土地利用/土地覆蓋遙感監(jiān)測數據分類系統(tǒng)Tab.1 Classification system of land use/land cover remote sensing monitoring data in China

2.2 MODHYDROLOG 模型簡介

MODHYDROLOG 模型是一個復雜的地下水補給模型，最初應用于澳大利亞［21，22］。該模型有5 個儲存庫和15 個參數。其結構及主要參數分別如圖2、表2 所示。MODHYDROLOG 模型模擬的徑流產出由3個部分組成：超滲地表徑流、蓄滿壤中流和基流。這3 種徑流在MODHYDROLOG 模型中的主要模擬過程為：首先降雨受到截留損失，超過最大截留儲量后剩余的降水一部分入滲進入土壤，另一部分形成超滲地表徑流和洼地蓄水。入滲的土壤水又會經過3個分配過程：補充地下水量，補充土壤含水量和蓄滿壤中流。以上過程中，截留量、洼地容量、土壤含水量都會受到蒸發(fā)作用，前兩部分蒸發(fā)按潛在蒸發(fā)能力計算，而在實際計算中，一般用實測蒸發(fā)量代替而作為模型的一個輸入量。模型原理及說明參見文獻［21，22］。

表2 MODHYDROLOG 模型參數及取值范圍Tab.2 Parameters and values range of MODHYDROLOG

圖2 MODHYDROLOG 模型結構Fig.2 Structure of MODHYDROLOG

2.3 參數優(yōu)選

上述模型用Nelder-Mead 算法進行參數優(yōu)化。Nelder-Mead算法是由Nelder和Mead于1962年在Spendley等［23］工作的基礎上設計出新的單純形搜索方法，該方法為求解無約束優(yōu)化問題的局部搜索算法，并無須目標函數的任何導數信息。對于D個變量的函數最小化問題，Nelder-Mead 單純形法使用反射、擴張、收縮、壓縮等操作，通過比較單純形的D+1 個頂點的目標函數值，用新的點取代目標函數值最大的頂點，逐步迭代并不斷更新，最終逼近問題的最優(yōu)解。本文在利用系列數據優(yōu)化評價的基礎上，對模型進行逐年優(yōu)化，推求不同年份的模型參數，以體現(xiàn)參數的時變特征。

2.4 模型評價指標

驗證模型準確性和評價模擬效果，可以通過以下統(tǒng)計學參數來反映。

主要有決定系數（R2），納什效率系數（NSE）［24］，標準偏差（Dv）［25］、克林效率系數（KGE）［26］：

式中：Qs和Qo分別代表第i 個樣本的模擬值和實測值；μs、μo分別為模擬值和實測值的平均值；σs、σo分別為模擬值和實測值的標準偏差；r為模擬值和實測值的相關系數。

兩種評價因子R2、NSE與KGE數值越接近1，則徑流的模擬值與實測值的差值越小，模擬效果越好。通常情況下，R2、NSE與KGE評價因子數值都在0.5 以上的徑流模擬就可表明模型精度較好［27］。據Donigian［28］的標準，認為Dv在觀測值的10%范圍內的模擬已經達到了“非常好”的標準，越接近于0 則表明模擬效果越好。

3 結果與分析

3.1 參數優(yōu)化及精度分析

率定期選擇漓江流域2008.1.1-2016.12.31 的降水、潛在蒸散發(fā)、氣溫等氣象數據，將2007全年數據年作為模型的預熱期。應用Nelder-Mead 算法對MODHYDROLOG 模型進行參數優(yōu)化選擇，設置待調參數上下邊界條件，以克林效率系數（KGE）為率定指標。模型優(yōu)化結果顯示，模型的多年KGE值接近于0.8，通過模擬效率系數檢驗；將最優(yōu)參數結果代入MODHYDROLOG模型，對2008-2016年漓江流域徑流過程進行模擬，比較分析決定系數（R2），納什效率系數（NSE），標準偏差（Dv），表明模型的決定系數R2、納什效率系數NSE均超過0.7，標準偏差Dv均在0.2%以內，達到了“非常好”的標準。因此，該模型對漓江流域日徑流模擬精度較高。圖3 給出了MODHYDROLOG 模型日徑流模擬結果，由圖3可知，兩種模型的日徑流過程與實測徑流過程具有很好的一致性，模型在枯水年、豐水年都能夠較好的再現(xiàn)日徑流變化過程，同時，洪峰的量級和出現(xiàn)時間基本一致，模擬效果令人滿意。

圖3 MODHYDROLOG 模型徑流模擬結果Fig.3 Runoff simulation results of MODHYDROLOG

3.2 蒸散發(fā)和徑流成分模擬結果

MODHYDROLOG 模型將蒸散發(fā)計算分為3層：植被截留蒸散發(fā)、土壤蒸散發(fā)和洼地蒸散發(fā)，總蒸發(fā)量為三層之和，計算得到的年蒸散發(fā)量見表3。

表3 MODHYDROLOG 模型模擬蒸散發(fā)情況 mmTab.3 Simulation of evapotranspiration by MODHYDROLOG

徑流成分是由水源劃分決定的，MODHYDROLOG 是三水源模型，由地表徑流、壤中流和基流3 種徑流成分構成。MODHYDROLOG 模型模擬所得各徑流量如圖4所示。2008-2016年地下水基準面、深層滲漏量、基流變化如圖5 所示，該圖體現(xiàn)了模擬期內地下水基準面年均值、深層滲漏量和基流年總值之間的相互關系。

圖4 模擬徑流成分比較Fig.4 Comparison of simulated runoff components

圖5 地下水基準面、深層滲漏、基流變化情況比較Fig.5 Comparison of changes of groundwater datum，deep leakage and base flow

3.3 土地利用變化與模型參數相關分析結果

將土地利用相關指標進行統(tǒng)計，通過綜合對比模型的相關參數來揭示徑流的變化特征，從而反映漓江流域的特征變化情況。如圖6 所示，將插值后的土地利用類型及變化與模型參數進行相關性分析。結果表明，不同土地利用類型的改變與參數的相關性差異較大。其中參數sub與草地為極顯著相關（p≤0.01）。參數dlev與城鄉(xiāng)，k1 與耕地和水域呈顯著相關（p≤0.05）。此外，參數crak與林地和草地、em與水域、dlev與耕地、林地、水域，k1與城鄉(xiāng)，k3與耕地、林地、城鄉(xiāng)相關（p≤0.1）。相關參數表明，土地利用變化通過影響不同徑流成分的比例影響產匯流過程，最終作用于模擬徑流量。主要是地下水儲量及其與河道之間的水量交換等相關參數與土地利用的變化密切相關，與植被作用相關的模型參數也有微弱相關性。而滲透（coeff、sq、vcond）和洼地的相關參數（dsc、ads、md）并無明顯顯著性。

圖6 土地利用變化與模型參數的相關性Fig.6 Correlation between land use change and model parameters

3.4 土地利用方式變化對流域降水-徑流特征的影響

3.4.1 耕地面積變化對流域降水-徑流特征的影響

圖7 表明，耕地面積與參數dlev 表現(xiàn)為顯著負相關，相關系數為-0.63，說明耕地會抑制地下水儲量基準面抬升，使得地下水深層滲漏量（與外流域的水量交換量）增多，不利于徑流的產生；同時，耕地面積與k1、k3 呈顯著正相關，相關系數分別為0.72 和0.59，其中與k1 為極顯著，說明耕地會促進地下水含水層與河流的水量交換作用，地下水補給河道水量，從而產生更多的徑流。此外，耕地面積與insc和em 分別呈微弱的正相關和負相關，反映出耕地面積增加有促進冠層截流和減緩植被蒸騰的作用，一定程度上增加土壤蓄水量和徑流量。

圖7 耕地與參數的相關性Fig.7 Correlation between cultivated land and parameters

3.4.2 林地面積變化對流域降水-徑流特征的影響

圖8 表明，林地面積與sub 呈微弱的正相關，因而林地對壤中流的形成起著微弱的促進作用，這是因為林地促進土壤非均質性增加更易形成壤中流［3］。林地面積與crak 呈顯著負相關，即表明林地會抑制入滲水量對地下水的優(yōu)先補給，即在土壤水分達到飽和前進入地下水的水量將會減少。林地面積對dlev、k1、k3 的影響與耕地相同，林地會降低地下水基準面，使得地下水深層滲漏增加，且更易于補給河道，有增加徑流的作用。

圖8 林地與參數的相關性Fig.8 Correlation between forest and parameters

3.4.3 草地面積變化對流域降水-徑流特征的影響

圖9 表明與耕地和林地不同的是，草地并沒有體現(xiàn)出與地下水相關參數的顯著性。唯一的是草地與crak 為正相關，相關系數為0.62，也就是說草地一定程度上促進了在水分入滲過程中補給地下水的比例。同時，草地與參數sub表現(xiàn)出極顯著負相關，相關系數為-0.84，即減小了壤中流的形成比例，這與增加地下水補給相一致，也體現(xiàn)了草地與林地的差別。另外，草地與土壤最大儲水量呈微弱正相關，顯著性系數為0.13，即一定程度上會促進土壤最大儲水量的增加。因此，草地較全面的體現(xiàn)了入滲水的3個分配方式的相關性，促進了土壤的儲水能力，增加了地下水的補給比例，降低了壤中流形成比例。

圖9 草地與參數的相關性Fig.9 Correlation between grassland and parameters

3.4.4 水域面積變化對流域降水-徑流特征的影響

圖10 表明水域面積與參數dlev呈顯著正相關，相關系數為0.59，說明水域面積會促進地下水基準面的抬升，減少地下水的深層滲漏量，有增加徑流的作用。與參數k1 和k3 呈顯著負相關，相關系數分別為-0.71 和-0.58，其中與k1 為極顯著相關性，顯著系數為0.03，表現(xiàn)出含水層基流交換減弱，并且補給河道的基流量減小，不利于徑流形成。此外，水域面積與參數insc表現(xiàn)出微弱的負相關，而與em呈顯著正相關，表明水域面積增加的同時會減小植被面積的比例，從而影響整個流域的截留能力，提供更多的水分供流域蒸發(fā)。

圖10 水域與參數的相關性Fig.10 Correlation between water areas and parameters

3.4.5 城鄉(xiāng)面積變化對流域降水-徑流特征的影響

圖11 表明城鄉(xiāng)面積與dlev呈極顯著正相關，相關系數為0.67，即城鄉(xiāng)面積有利于抬升地下水的基準面，減少地下水的深層滲漏量，增加河道徑流量。城鄉(xiāng)面積與參數k1、k3 呈顯著負相關，表現(xiàn)出弱化含水層基流對河道的補給，不利于徑流形成。城鄉(xiāng)面積與insc呈反比，說明對冠層截留有抑制作用，可適當增加凈雨，進而形成較多的徑流量。

圖11 城鄉(xiāng)面積與參數的相關性Fig.11 Correlation between urban and rural areas and parameters

總體來說，5 種土地利用類型中除草地外，其余4 種土地利用類型方式都對地下水補給施加不同程度的影響，進而影響基流的水量交換，通過參數dlev、k1、k3 從深層滲漏和地下水補給河流得以表現(xiàn)。土壤滲透的相關參數（coeff、sq、vcond）以及洼地的相關參數（dsc、ads、md）與流域降水-徑流特征相關性不顯著，即土地利用方式的改變并不會顯著地通過改變土壤滲透性和流域內洼地的面積來影響徑流形成。因此，土地利用類型對不同水文過程的影響程度不同，因而單一組分對產匯流的影響也存在差異，表4歸納了土地利用類型對相關參數及徑流的影響。

表4 土地利用類型對相關參數及徑流的影響Fig.4 Effects of land use types on relevant parameters and runoff

4 討 論

研究結果表明，模型參數經率定后，MODHYDROLOG 模型能較好地模擬漓江流域的日徑流過程，NSE 和R2均超過0.7，Dv在0.2%以內。主要由于MODHYDROLOG 模型考慮了超滲產流和蓄滿產流兩種產流機制，具有一定的物理基礎，在枯水期和豐水期都能很好的再現(xiàn)日徑流的變化過程，因而能在處于濕潤地區(qū)的漓江流域適用性較好，一定程度上反映了研究區(qū)漓江流域的產匯流特征。但由于漓江流域全年降雨量較大且集中降雨較多，這使得模型主要以超滲產流的機制為主，因此在枯水期存在一定誤差。

當耕地、林地面積的減小，草地、水域面積的增加時，與之對應的是漓江流域的地下水基準面的抬升，深層滲漏減少，同時基流交換愈加強烈，地下水大量補充河道產生徑流。如圖5所示，由于地下水基準面抬升作用，深層滲漏量有減小趨勢；同時，地下水與河道間的基流交換有增強趨勢，產生更多徑流。林地的土壤類型多為不飽和雛形土和疏松巖性土，土質較為疏松，經土壤入滲補充地下水的能力較強，所以地下水補給量高，具有較強的水土保持和涵養(yǎng)水源的功能，也因根系作用易于形成壤中流。草地的土壤類型多為飽和黏性土，此類土表層疏松，但深層的土層黏粒含量相對較高，這使降水通過入滲來補給地下水變得困難，而且草地的根系分布較淺，僅能從淺層土壤中吸收水分，因而與地下水有關參數的相關性并不顯著。研究時段內草地面積的變化僅微弱上升并不顯著，因此并不能作為反映流域特征變化的依據［29］。另外，水域面積的增加造成更多的水量直接匯入河道，同時加大了下滲后介于地下水和河道之間的基流交換；與之同樣效果的是，城鄉(xiāng)建設增加了不透水面積，因此更多的降水直接進入河道形成徑流。

模擬結果表明，模型中蒸散發(fā)和壤中流等相關參數與大部分土地利用指標相關性較弱，一方面源于植被的作用，另一方面也體現(xiàn)了土地利用方式的改變對漓江流域的蒸散發(fā)和壤中流等影響不大。此外，土地利用指標與模型參數的相關性分析結果表明，土地利用方式與滲透和洼地的相關參數（coeff、sq、vcond、dsc、ads、md）的顯著性水平較低（p＞0.2），并無表現(xiàn)出明顯的相關性，與Francis Chiew 認為這類參數與流域特征相關性極弱的研究結果相一致［3］。

通過土地利用變化與模型參數的相關分析，針對不同土地利用類型作用于模型的時變參數的情況，最終解釋土地利用變化對流域產匯流過程的影響，即流域降水-徑流特征。土地利用變化對流域徑流量的影響取決于流域內土地利用變化的程度，由于本文所采用的時間序列較短，土地利用指標的精度不高，可能存在參數與流域特征指標顯著性不高的情形。另外，對于參數如k1、k3 等雖與土地利用指標存在顯著性，但由于其作為方程的系數且物理意義不太明確，故在解釋流域的具體特征變化上略顯無力。此外，模型對于地表水與地下水的相互作用能較好的概括，但不一定對流域特征較好地表達；流域中有新建水庫參與降水徑流過程，但參數中體現(xiàn)不足，可能與運行時間較短有關；在分析土地利用方式變化時，沒有考慮不同方式之間的轉化，僅單一地考慮某一利用土地面積的變化造成的降水-徑流特征的影響。

模型中與地下水相關的參數與土地利用類型指標變化有較強的相關性，而這些參數又將直接作用于降雨轉化為徑流的整個過程，因而一定程度上解釋了漓江流域降雨-徑流特征。具體而言，土地利用類型的變化主要與模型中補給地下水、地下水基準面及基流（crak、dlev、k1、k3）有較強的相關性。在概念性水文模型中，如dlev 一類的參數，不僅反映地下水的深層滲漏，也對模擬流域水量平衡方面發(fā)揮補償調節(jié)作用［30］，即有些土地利用改變，會造成流域降水徑流平衡誤差增大時，模型以dlev 參數來進行調節(jié)補償。另一方面，研究區(qū)中下游平原區(qū)多石灰?guī)r分布，因此可以通過地下暗河發(fā)生地下水深層滲漏。

5 結 論

（1）MODHYDROLOG 模型能很好的模擬漓江流域的徑流量變化情況，但由于漓江流域降雨充沛且分配不均，主要以超滲產流為主、蓄滿產流為輔，同時在降雨強度較小的枯水期容易高估徑流量，因而產生偏差。

（2）相關分析表明，模型參數主要是地下水相關參數dlev、k1、k3 與絕大部分土地利用類型變化關系密切，個別土地類型的變化與植被相關參數sub、crak、em 有著不同程度的相關性。主要表現(xiàn)為耕地和林地面積與參數dlev 負相關，與k1、k3 正相關，水域和城鄉(xiāng)反之。此外，林地和草地與參數crak 相關，水域與參數em 也達到了相關性。而滲透（coeff、sq、vcond）和洼地的相關參數（dsc、ads、md）并無明顯顯著性。

（3）對模型參數和下墊面指標進行相關性分析，分析土地利用對產匯流徑流成分的影響，主要體現(xiàn)在耕地和林地抑制地下水儲量基準面抬升，使得地下水深層滲漏量增多，同時促進了含水層與河道之間的基流交換作用，有利于地下水對河道的補給。草地可促進水分入滲補給地下水的比例，減小了壤中流的形成比例，促進土壤的儲水能力。水域和城鄉(xiāng)與耕地和林地的作用機制相反，其對基流交換的抑制作用減緩了地下水對河道的補給量。