新安江模型在資料匱乏的長江中下游山區(qū)中小流域洪水預(yù)報應(yīng)用*

龔珺夫,陳紅兵，朱 芳，李永凱，崔 明，李致家

(1：河海大學(xué)水文水資源學(xué)院，南京 210098)(2：宜昌市水文水資源勘測局，宜昌 443000)

可靠的洪水預(yù)報是防洪決策中的重要依據(jù)[1]，是減少甚至避免洪水破壞的第一道防線. 在水文學(xué)研究及實際工程應(yīng)用中，相較于水文資料充足的地區(qū)，更常見且更有挑戰(zhàn)性的是對資料匱乏流域進行洪水預(yù)報. 國際水文學(xué)會(IAHS)曾在20世紀初提出了水文資料匱乏流域的洪水預(yù)報(PUBs)10年計劃[2]，期間涌現(xiàn)了大量研究成果[3-5]. 資料匱乏流域的洪水預(yù)報研究一般被概括為以下兩部分[6]：(1)直接針對徑流或徑流指標(biāo)進行區(qū)域化；(2)針對水文模型參數(shù)進行區(qū)域化. 兩者都是基于回歸方法或目標(biāo)流域與測量流域之間的某種距離度量. 首先，簡單的回歸模型顯示出較高的應(yīng)用價值，如Mazvimavi等[7-8]在非洲52個流域建立了徑流指標(biāo)與流域?qū)傩缘亩嘣€性回歸模型，以期通過年平均降雨、流域平均坡度、河網(wǎng)密度和土地覆蓋等流域?qū)傩怨烙嫙o資料流域的年平均徑流量、基流指數(shù)、歷時曲線和月平均徑流量等徑流指標(biāo). Maréchal等[9]基于英國土壤類型水文計劃(HOST)的土壤數(shù)據(jù)，使用了一個簡單的回歸方法較精確地預(yù)測了未測量流域的基流指數(shù)等徑流指標(biāo). 也有部分研究基于更復(fù)雜的回歸模型，如Laaha等[10-11]揭示了根據(jù)區(qū)域特點定制回歸模型而不是使用全局模型的好處，其研究表明基于按季節(jié)性分組的流域回歸模型提供了魯棒性更強的區(qū)域化結(jié)果. 同樣地，地質(zhì)統(tǒng)計學(xué)方法也被用于估計資料匱乏流域的水文變量，比如Sk?ien等[12]提出了拓撲克里金法，可以用于解釋水動力及地貌特征的空間離散性，該方法避免了輸入數(shù)據(jù)的誤差和參數(shù)可識別性問題. 部分研究[13]認為該方法提供了比回歸模型更穩(wěn)健的估計. 另一方面，研究者們越來越認識到，流域空間上的接近不一定導(dǎo)致流域產(chǎn)流模式的相似[14]，當(dāng)使用更具有水文學(xué)意義上的距離量度時，基于距離的方法的估計效果會顯著提高[15]. 如Merz等[16]將拓撲克里金法與流域特征相結(jié)合，提高了該方法的預(yù)測性能. Dornes等[17]認為只要保證流域之間的水文相似性，水文模型的參數(shù)甚至可以傳遞數(shù)千千米. 時至今日，資料匱乏地區(qū)的水文預(yù)報獲得了長足進步，但仍有大量難題無法解決[18].

我國有防洪任務(wù)的重點中小河流(匯流面積小于3000 km2)有5186條[19]. 近年來，受氣候變化影響，由局地強降水造成的中小河流突發(fā)性洪水頻繁發(fā)生，已成為造成人員傷亡的主要災(zāi)種. 據(jù)統(tǒng)計，我國中小流域澇災(zāi)害和山洪地質(zhì)災(zāi)害損失約占全國洪澇災(zāi)害經(jīng)濟損失的70%～80%，對人民群眾的生命財產(chǎn)安全構(gòu)成了嚴重的威脅[20]. 由于大部分中小河流源短流急，洪水具有歷時短、上漲快的特點，加之水文資料短缺，傳統(tǒng)的預(yù)報方案通常難以對這些地區(qū)的洪水做出有效的預(yù)報[21]. 區(qū)域化方法是解決此問題的思路之一，目前針對無資料中小流域的區(qū)域化研究較少，姚成等[22]使用區(qū)域化方法(空間臨近法和回歸法)將API 模型和新安江模型應(yīng)用于大別山區(qū)及皖南山區(qū)的29個中小流域，得到了較高的預(yù)報精度. 但是在現(xiàn)有的區(qū)域化研究中，測量流域與目標(biāo)流域大多空間分布均勻且空間距離較近. 在這些地區(qū)，基于空間距離的區(qū)域化方法(以下稱為空間臨近法)由于操作簡單且所需資料少而多被使用. 同時，基于水文屬性距離的區(qū)域化方法(以下稱為物理相似法)和回歸方法也被用來提高預(yù)測性能. 在自然條件與歷史因素的限制下[23]，設(shè)立在我國濕潤山區(qū)中小流域的水文測站分布稀疏，可能導(dǎo)致某一具體區(qū)域內(nèi)的測量流域數(shù)量較少. 是否可以通過其他資料充足但相隔較遠地區(qū)的測量流域作為參證流域進行參數(shù)估計？為了回答此問題，本研究關(guān)注于目標(biāo)流域與測量流域空間位置較遠時，區(qū)域化方法是否可以有效估計參數(shù). 此時，測量流域與目標(biāo)流域之間顯著的空間差異導(dǎo)致空間臨近法無法使用，物理相似法及回歸方法是否適用，需要接下來進一步的研究. 本研究將皖南山區(qū)與鄂西山區(qū)分別作為已觀測區(qū)及未觀測區(qū)，選用新安江模型進行洪水模擬，并驗證區(qū)域化方法(物理相似法及回歸法)在研究區(qū)的適用性. 正如Prieto等[24]所言，單獨使用參數(shù)區(qū)域化有諸多弊端，本研究基于新安江模型參數(shù)的物理意義先估算出部分參數(shù)，剩余參數(shù)使用區(qū)域化方法(物理相似法及回歸法)估計，并比較不同參數(shù)估計方案的效果，提出在研究區(qū)內(nèi)最為適用的參數(shù)估計方案，為長江中下游資料匱乏的山區(qū)中小流域洪水預(yù)報提供參考.

1 研究區(qū)概況及資料來源

研究區(qū)域為長江中下游地區(qū)的32個山區(qū)中小流域，其中位于皖南山區(qū)的29個流域為有觀測資料的測量流域，鄂西山區(qū)的3個流域被視為資料匱乏的目標(biāo)流域，研究區(qū)域見圖1. 各流域面積均不超過1000 km2，均為自然閉合流域，受水工建筑物的影響較小，流域平均海拔多在400～1000 m，其中霧渡河流域平均海拔最高，為1109 m. 各流域年平均降水量均在1200～1500 mm，年平均蒸發(fā)量均在600～800 mm，年平均氣溫14～16℃. 研究流域大多地勢陡峭，河網(wǎng)密集，平均坡度大多在10～30°，河網(wǎng)密度大多在0.6～1.2 km/km2. 研究區(qū)內(nèi)土壤類型以壤土和黏土為主，土壤飽和滲透系數(shù)大多在3～12 mm/h. 土地利用(LUCC)以林地為主，大部分流域林地面積占比70%以上，甚至部分達到100%. 研究區(qū)內(nèi)植被類型以木本植物為主，大多為常綠針葉林、常綠闊葉林及灌木.

圖1 研究區(qū)地理位置

降雨、徑流及蒸發(fā)資料來自當(dāng)?shù)厮木郑瑫r間步長均為1 h. 29個測量流域均具有至少30年的觀測資料. 3個目標(biāo)流域中，霧渡河有2008-2017年間共11場觀測洪水，茅坪河有2016-2018年間共5場觀測洪水，下牢溪僅有2016年共2場觀測洪水. DEM來自中國科學(xué)院地理空間數(shù)據(jù)云，分辨率為30 m，用于提取流域地形、地貌屬性；土地利用數(shù)據(jù)來自中國科學(xué)院遙感解譯的2010年LUCC圖，分辨率為1 km，用于提取流域土地利用屬性；土壤數(shù)據(jù)來自于HWSD(世界土壤數(shù)據(jù)庫)的土壤類型圖，分辨率為1 km，用于提取流域土壤屬性.

2 研究方法

本研究的研究思路大致如下：首先在測量流域進行新安江模型的校準，得到測量流域的新安江模型校準參數(shù)集. 然后進行流域物理屬性的篩選，依據(jù)篩選出的流域物理屬性集，運用區(qū)域化方法(回歸法、物理相似法)，將測量流域的校準參數(shù)轉(zhuǎn)移至目標(biāo)流域. 與此同時，使用基于參數(shù)物理意義的估算方法(以下簡稱物理估算法)，估算出部分參數(shù). 最后，將物理估算法與區(qū)域化方法相結(jié)合，提出針對目標(biāo)流域的多個綜合參數(shù)估計方案，并比較不同參數(shù)估計方案的效果，得到在研究區(qū)內(nèi)最為適用的參數(shù)估計方案. 本研究的技術(shù)路線圖見附圖Ⅰ.

2.1 新安江模型及物理估算法

2.1.1 新安江模型 新安江模型是趙人俊[25]于1984年前后提出的概念性模型，在我國濕潤及半濕潤地區(qū)具有較高的預(yù)報精度，已得到廣泛應(yīng)用[26-28]. 該模型共有16個參數(shù)，各參數(shù)物理意義見文獻[29]，其中不敏感參數(shù)10個(WM、WUM、WLM、B、C、IMP、EX、CI、KE、XE)，敏感參數(shù)7個(KC、SM、KI(KG)、CG、CS、L). 根據(jù)趙人俊等[30]提出的客觀優(yōu)選法，不敏感參數(shù)不參與率定，直接賦予固定值. 敏感參數(shù)的率定使用人工優(yōu)選結(jié)合SCE-UA算法[31]，這樣既可以快速穩(wěn)定的優(yōu)選參數(shù)，又能有效避免自動優(yōu)選參數(shù)可能不合理的問題. 選用洪峰相對誤差、峰現(xiàn)時間誤差，納什效率系數(shù)(NSE)對參數(shù)率定結(jié)果進行評價. 依據(jù)《水文情報預(yù)報規(guī)范》[32]，洪峰相對誤差在20%以內(nèi)，峰現(xiàn)時間誤差在3 h以內(nèi)為合格.

在目標(biāo)流域，不敏感參數(shù)同樣直接賦予固定值，敏感參數(shù)需要進行估計. 在敏感參數(shù)中KI與KG滿足線性約束，只需估計KI即可，因此本研究在目標(biāo)流域需要估計的參數(shù)為6個，見表1. 在新安江模型的待估計參數(shù)中，壤中流出流系數(shù)KI和河網(wǎng)水消退系數(shù)CS可以基于參數(shù)的物理意義進行估算(物理估算法).

表1 目標(biāo)流域待估計參數(shù)

2.1.2 壤中流與地下水出流系數(shù)KI與KG土壤質(zhì)地可以用來推斷土壤的性質(zhì)，進而估算自由蓄水量對壤中流及地下水的出流系數(shù)(新安江模型中的KI與KG).KI和KG之和決定了自由水的排水速率，趙人俊等[29]認為KI+KG的值可固定為0.7；KI/KG決定了壤中流與地下水出流的比例，姚成等[33]在GXM模型中提出了針對柵格的KI/KG估計公式，將其修改成適用于新安江模型的形式：

(1)

表2 不同質(zhì)地土壤凋萎含水量

式中，mr為自由水出流校正系數(shù)，根據(jù)前人的研究成果[33-35]，在濕潤流域mr可取為1；θwp為土壤凋萎含水量，以體積含水量表示，可根據(jù)土壤質(zhì)地取值，見表2[36].

2.1.3 河網(wǎng)水消退系數(shù)CS李致家等[37]針對二叉樹結(jié)構(gòu)河網(wǎng)，對河鏈進行編號，鏈接河源的河鏈編號記為r=1，下游河鏈的編號等于該條河鏈的上游總鏈數(shù)+1. 并根據(jù)矩形河道線性蓄泄關(guān)系(式(2))，引入一個無量綱時間tα：

(2)

(3)

式中，Q(r,t)為河鏈r的出流量，m3/s；W(r,t)為河鏈r的蓄量，m3；v為平均河道流速，m/s；L為平均河鏈長，m.

在資料匱乏地區(qū)一般沒有實測流速資料，可以使用經(jīng)驗公式估算. Rodríguez-Iturbe等[38]基于運動波理論提出如下流速公式：

(4)

(5)

式中，p為次洪的平均雨強，cm/h；A為流域面積，km2；B為平均河寬，m；S0為平均河道坡度；n為河道曼寧糙率系數(shù)，根據(jù)已有研究[39]，這里可取n=0.025.

李致家等[37]在此基礎(chǔ)上推導(dǎo)出遞推形式的河鏈蓄量公式：

(6)

式中，r1與r2分別為河鏈r上游的兩條河鏈.

河網(wǎng)消退系數(shù)CS表示河網(wǎng)退水速度，可以用河網(wǎng)出口處計算時段起止的流量比表示[40]：

(7)

式中，Qr(t)為河網(wǎng)出口流量，m3/s；ts表示計算時段起始時刻；te表示計算時段終止時刻.

將式(2)、(3)代入式(7)，得到CS與蓄量的關(guān)系[41]

(8)

為了避免使用單一時段造成的估計誤差，通過計算N1個時段的起止蓄量，使用最小二乘法估計CS的值[41]：

(9)

2.2 參數(shù)區(qū)域化方法

區(qū)域化是在水文資料匱乏流域進行水文預(yù)報的有效方法，所謂區(qū)域化即是將水文模型參數(shù)值從水文資料豐富的測量流域轉(zhuǎn)移到水文資料匱乏的目標(biāo)流域的過程[42]. 目前，常用的區(qū)域化方法主要有物理相似法、回歸法和空間鄰近法. 物理相似法[43]將參數(shù)值從測量流域轉(zhuǎn)移到目標(biāo)流域，該測量流域的屬性(氣候與下墊面等)與目標(biāo)流域的屬性相似. 回歸法[44]則是在測量流域上建立校準參數(shù)與流域?qū)傩?氣候與下墊面等)之間的關(guān)系，然后根據(jù)目標(biāo)流域的流域?qū)傩院鸵呀⒌幕貧w關(guān)系估計目標(biāo)流域的參數(shù)值. 空間鄰近法[45]假定空間上鄰近的流域，下墊面與氣候條件也應(yīng)該相似，因此目標(biāo)流域的參數(shù)值通常來自一個或多個附近的測量流域.

2.2.1 物理相似法 單流域物理相似法選擇與目標(biāo)流域水文物理相似性最高的測量流域作為參數(shù)轉(zhuǎn)移的供體(以下簡稱參證流域)，將參證流域的參數(shù)值直接移用到目標(biāo)流域. 利用秩累積相似性[46]作為流域物理相似性指標(biāo)，對于流域?qū)傩约?區(qū)域化研究中被選擇的流域?qū)傩缘募?里的每個屬性，具有與目標(biāo)流域最相似屬性的測量流域的秩記為1，具有第二相似屬性的測量流域的秩記為2，依此類推. 計算流域?qū)傩约锩總€屬性的秩號，選擇總秩最小的測量流域作為參證流域，每個用于區(qū)域化的屬性在秩累積中都具有同等的權(quán)重. 物理相似法要討論的另一個重要問題是參證流域的數(shù)量，傳統(tǒng)物理相似法僅考慮一個最相似的流域，但正如Mcintyre等[47]指出的，來自最相似流域的參數(shù)集不一定產(chǎn)生最佳結(jié)果. Viviroli等[48]建議選取前5個相似流域作為參證流域，可以使計算效率與模型準確性達到最佳平衡. 因此本研究提出多流域物理相似法，選取秩累積最小的前n個測量流域作為參證流域，將每個參證流域的參數(shù)集分別移用至目標(biāo)流域，然后在每個計算時間步長上取模擬流量的中位數(shù)，組成最終的模擬結(jié)果. 使用中位數(shù)是為了避免模擬結(jié)果出現(xiàn)過度平滑現(xiàn)象.

本研究為了定量分析物理相似法的適用性，引入用戶加權(quán)的歐氏距離[48]作為相似性度量，用戶加權(quán)的歐氏距離越小，說明兩個流域間的水文物理相似性越高：

(10)

式中，Dφ(x,y)為x、y兩流域間的用戶加權(quán)的歐氏距離；N2為流域?qū)傩约镌貍€數(shù)；ATTα(x)為流域x的第α個標(biāo)準化流域?qū)傩?；φα為用戶賦予第α個標(biāo)準化流域?qū)傩缘臋?quán)重，這里取φα=1[48].

2.2.2 回歸法 回歸法適用于測量流域數(shù)量較多(大于10個)且流域水文物理特征相似的地區(qū)[49]. 回歸法假設(shè)流域特征與模型參數(shù)間存在較好的相關(guān)關(guān)系，因此具有一定的局限性[50]. 部分模型參數(shù)可能與流域?qū)傩韵嚓P(guān)性較低，對于這類參數(shù)，難以建立可靠的回歸關(guān)系，不適用回歸法. 另外，回歸法中流域?qū)傩缘倪x擇十分重要，若流域?qū)傩赃x擇不恰當(dāng)，回歸法的預(yù)測性能較差. 在以往的研究中[51-53]，大多根據(jù)專家經(jīng)驗人工選擇參與回歸的流域?qū)傩?，具有極大的主觀性；同時為了方便處理，通常使用線性回歸，而線性回歸很難代表實際水文情況[54]，具有較大缺陷. 本研究引入Boruta算法[55]進行流域?qū)傩院Y選，以減少流域?qū)傩赃x擇的主觀性，并使用隨機森林算法[56]進行非線性回歸.

隨機森林是bagging算法[57]的一個變體，以決策樹為基學(xué)習(xí)器，構(gòu)建多顆決策樹來組成一個“森林”. 隨機森林使用自助采樣法(bootstrap sampling)[58]，以包含β個樣本的數(shù)據(jù)集構(gòu)建決策樹的時候，依次有放回的取出一個數(shù)據(jù)放入采樣集，重復(fù)取樣β次，這樣在創(chuàng)建多顆決策樹時就有部分數(shù)據(jù)被多次選擇，同時也有部分數(shù)據(jù)始終沒有被選擇，未被選擇的數(shù)據(jù)大約占數(shù)據(jù)總量的36.8%[56]，這部分數(shù)據(jù)被稱為袋外數(shù)據(jù)(out-of-bag data). 同時，隨機森林隨機地從θ個自變量中選擇部分自變量進行決策樹節(jié)點的確定，因此每次構(gòu)建的決策樹都可能不一樣. 最終隨機森林預(yù)測結(jié)果為“森林”中每決策樹單獨預(yù)測結(jié)果的期望值. 隨機森林相比于其它回歸方法有以下優(yōu)點：(1)可以處理變量的多元共線性問題；(2)模型結(jié)構(gòu)十分靈活，不需要預(yù)先確定回歸方程的形式；(3)隨機性的引入使得隨機森林不容易過擬合，有很好的抗噪聲能力；(4)由于存在袋外數(shù)據(jù)，使隨機森林沒有必要進行交叉驗證或單獨留出驗證集，袋外誤差(out-of-bag error)就是模型泛化誤差的一個無偏估計[59]，其結(jié)果近似于k折交叉驗證.

Boruta算法是一個基于隨機森林的特征選擇方法，使用Z分數(shù)作為特征重要性度量，考慮了隨機森林中決策樹之間平均準確度損失的波動. 該方法將每個原始變量進行隨機復(fù)制，創(chuàng)造出原始變量的“陰影”變量，將“陰影”變量與原始變量合并，并使用隨機森林對這個擴展的信息系統(tǒng)分類，收集并計算出Z分數(shù). 當(dāng)一個真實變量比最好的“陰影”變量具有更高的Z分數(shù)，則認為其對回歸模型是有意義的. 每次迭代中，不斷刪除不重要的特征，直到所有真實變量得到確認或拒絕. Boruta算法的具體原理可見Kursa等[60]，此處不再贅述.

3 流域?qū)傩院Y選及方法適用性評價

3.1 流域?qū)傩院Y選

選取適當(dāng)?shù)牧饔驅(qū)傩允沁M行區(qū)域化研究的必要前提，目前已有大量流域?qū)傩员惶岢鯷3,15,61-65]，這些研究普遍認為流域?qū)傩员仨氁罁?jù)不同的研究方法和研究區(qū)域進行選擇，不應(yīng)直接套用. 流域?qū)傩缘倪x擇不宜過多[66]，如英國洪水估算手冊[67]推薦使用3種流域?qū)傩詡鬟f水文模型參數(shù)；Samaniego等[68]在研究中使用了7種流域?qū)傩? 本研究依據(jù)研究區(qū)可獲得的數(shù)據(jù)情況、新安江模型參數(shù)的物理意義及山區(qū)中小河流的特點，從氣象、地形地貌特征、土地利用及土壤等方面選取25種流域?qū)傩?，見?. 首先選取流域基本屬性及流域的形狀屬性，包括流域面積、高程，長度、寬度及形狀系數(shù)等，這些屬性可能影響蒸散發(fā)折算系數(shù)KC等模型參數(shù). 其次，研究流域大多地勢陡峭、河網(wǎng)密集，反映流域坡度、河道坡度及河網(wǎng)密度的流域?qū)傩?如流域平均坡度、河道平均坡度、河網(wǎng)密度等)應(yīng)該被考慮，且這些屬性對流域匯流影響較大，可能影響河網(wǎng)水流消退系數(shù)CS、地下水消退系數(shù)CG和河網(wǎng)匯流滯時L等模型參數(shù). 地形指數(shù)是Beven等[69]提出的一個地形特征參數(shù)化指標(biāo)，是流域徑流源面積和地下水水位空間分布特征的近似表征[70]，本研究將其作為一個水文綜合指標(biāo). 土地利用和土壤類型綜合反映流域下墊面的產(chǎn)匯流條件，影響著自由水容量SM及壤中流出流系數(shù)KI等模型參數(shù).

通過皮爾遜相關(guān)分析剔除顯著強相關(guān)的流域?qū)傩?，篩選出最不相關(guān)的流域?qū)傩越M成流域?qū)傩约?，其中各個屬性的皮爾遜相關(guān)系數(shù)P<0.7且顯著性水平α>0.05[71]. Viviroli等[48,72]研究認為，在篩選出的流域?qū)傩约?，分別選取與每個待估計水文模型參數(shù)相關(guān)性最高的兩個流域?qū)傩?，組成最終的流域?qū)傩约?，區(qū)域化結(jié)果最優(yōu). 新安江模型的待估計參數(shù)有6個(KC、SM、KI、CG、CS、L)，通過皮爾遜相關(guān)分析，分別選取與待估參數(shù)相關(guān)性最高的兩個流域?qū)傩?，組成最終本研究使用的流域?qū)傩约?表3中帶*屬性). 流域?qū)傩约锇ㄒ韵?種流域?qū)傩?流域?qū)傩圆蛔?2個是因為有部分流域?qū)傩酝瑫r與多個模型參數(shù)高度相關(guān))：流域面積AREA、流域平均高程ELEVATION、河道平均坡度CSLOPE、面積坡度ASLOPE、河網(wǎng)密度DD、林地面積比FST和土壤飽和滲透系數(shù)SOLK.

表3 流域?qū)傩?/p>

圖2 流域物理相似性分析

3.2 適用性評價

本研究選取的模擬時間步長為1 h. 29個測量流域的新安江模型模擬結(jié)果可見文獻[22]，平均洪峰合格率、平均峰現(xiàn)時間合格率和平均確定性系數(shù)分別為84%、91%和0.84，表明在本研究區(qū)域，新安江模型的次洪模擬精度較高，是適用的水文模型.

測量流域與目標(biāo)流域空間距離較遠，主要分為皖南和鄂西兩大簇，因此空間近似法不適用. 但兩者具有相似的氣象條件，且研究流域均是山區(qū)中小河流流域，具有相似的下墊面條件，為使用物理相似法及回歸法提供了條件. 根據(jù)公式(10)，在流域?qū)傩缘膎維空間中，依次選取一個測量流域作為參證流域，分別計算與剩余的28個測量流域之間的歐氏距離，及與3個目標(biāo)流域之間的歐氏距離，比較結(jié)果見圖2. 可見兩者的均值相近，且目標(biāo)流域與測量流域之間的歐氏距離最大值未超過測量流域之間歐式距離的最大值，因此認為目標(biāo)流域與測量流域之間物理相似性較高，使用物理相似法是合理的.

對于回歸法而言，使用Boruta算法分別計算流域?qū)傩詫γ總€敏感參數(shù)的重要性，以評估其通過隨機森林建立回歸關(guān)系的可靠性. 結(jié)果顯示，只有河網(wǎng)水消退系數(shù)CS及河網(wǎng)匯流滯時L與流域?qū)傩约锏牟糠謱傩杂休^密切的關(guān)系，見圖3. 圖中藍色箱體代表“陰影”屬性Z分數(shù)的最小值、均值與最大值，大于“陰影”屬性Z分數(shù)最大值的流域?qū)傩员粯?biāo)記為“接受”，即綠色箱體，剩余的流域?qū)傩员粯?biāo)記為“拒絕”，即紅色箱體. 對剩余敏感參數(shù)而言，流域?qū)傩约锏乃袑傩缘腪分數(shù)均小于“陰影”屬性最大值，被Boruta算法標(biāo)記為“拒絕”. 即流域?qū)傩耘c這些參數(shù)沒有明顯相關(guān)性，相當(dāng)于隨機變量，因此無法使用回歸法對這部分參數(shù)進行區(qū)域化.

圖3 部分新安江模型參數(shù)的Boruta分析

CS與L分別與被接受的流域?qū)傩越⒒貧w關(guān)系. 在泛化誤差分析中，隨機森林使用袋外誤差作為泛化誤差，Wolpert 等[59]已經(jīng)證明隨機森林的袋外誤差結(jié)果近似于k折交叉驗證. 線性回歸及逐步回歸模型使用留一法(leave-one-out cross validation，即n折交叉驗證，n為樣本數(shù))估計泛化誤差. 以平均絕對誤差(MAD)、均方誤差(MSE)及偽決定系數(shù)(RSQ)[73]作為誤差評價指標(biāo)，其中RSQ取值在0～1之間，表示自變量對因變量的解釋程度，越接近1表示解釋程度越高：

(11)

式中，var(y)為因變量實際值的方差.

表4統(tǒng)計了將流域?qū)傩约锶苛饔驅(qū)傩宰鳛樽宰兞康木€性回歸模型、姚成等[22]在相同區(qū)域建立的逐步回歸模型和本研究建立的隨機森林回歸模型的泛化誤差. 隨機森林模型中CS的MAD、MSE與線性回歸差別不大，但是隨機森林模型中L的MAD、MSE明顯低于線性回歸模型. 另一方面，隨機森林模型中CS及L的MAD、MSE均小于逐步回歸. 說明在本研究區(qū)，隨機森林泛化能力優(yōu)秀，是較為適用的方法. 且隨機森林法的RSQ分別為0.46和0.58，均高于對應(yīng)的線性回歸方法和逐步回歸方法，表明Boruta方法篩選的流域?qū)傩耘cCS及L關(guān)系密切，建立的回歸模型較為可靠.

表4 回歸模型泛化誤差分析

4 方案設(shè)計與結(jié)果

設(shè)計了10種不同的方案進行目標(biāo)流域的參數(shù)估計：(1) 所有敏感參數(shù)取29個測量流域的平均值：(2) 所有敏感參數(shù)使用單流域物理相似法估計；(3) 所有敏感參數(shù)使用多流域物理相似法(2個參證流域)估計；(4) 所有敏感參數(shù)使用多流域物理相似法(5個參證流域)估計；(5)CS及L使用回歸方法，剩余參數(shù)使用單流域物理相似法估計；(6)CS及L使用回歸方法，剩余參數(shù)使用多流域物理相似法(2個參證流域)估計；(7)CS及L使用回歸方法，剩余參數(shù)使用多流域物理相似法(5個參證流域)估計；(8)L使用回歸方法估計，CS及KI(KG)使用物理估算法，剩余參數(shù)使用單流域物理相似法估計；(9)L使用回歸方法估計，CS及KI(KG)使用物理估算法，剩余參數(shù)使用多流域物理相似法(2個參證流域)估計；(10)L使用回歸方法估計，CS及KI(KG)使用物理估算法，剩余參數(shù)使用多流域物理相似法(5個參證流域)估計. 方案1是不進行任何優(yōu)化的基準參照組；方案2～4單獨使用了物理相似法；方案5～7結(jié)合了物理相似法和回歸法；方案8～10結(jié)合了物理相似法、回歸法和物理估算法. 方案2～4之間的區(qū)別是物理相似法中選取的參證流域數(shù)目分別為1、2、5個，方案5～7、方案8～10同理. 如此設(shè)計方案的目的是：(1) 比較不同方法或不同方法組合對參數(shù)估計精度的影響；(2) 比較不同參證流域數(shù)量對參數(shù)估計精度的影響.

使用平均洪量相對誤差、平均洪峰相對誤差(為了避免正負誤差抵消，均以絕對值計算)和納什系數(shù)(NSE)作為新安江模型效率的評價指標(biāo)，并將不同方案估計出的參數(shù)與率定的參數(shù)進行比較，各指標(biāo)的變化即是參數(shù)估計方案導(dǎo)致的模型效率降低：

(12)

表5統(tǒng)計了各方案估計參數(shù)的精度. 方案1作為不進行任何優(yōu)化的基準參照組，參數(shù)估計效果最差，在3個目標(biāo)流域NSE均遠小于0，且模型效率損失最大，說明使用此方案估計參數(shù)的新安江模型不可信. 方案2～10使用不同的方法進行了優(yōu)化，模型效率損失均比方案1有所降低，說明這些方案不同程度的提高了資料匱乏流域的參數(shù)估計精度. 比較三組選擇不同參證流域數(shù)量的物理相似法(方案2、3、4；方案5、6、7及方案8、9、10)，在茅坪河與霧渡河，模型效率損失為方案2<方案3<方案4，方案5<方案6<方案7，方案8<方案9<方案10. 在下牢溪，除方案8模型效率損失顯著較低外，其他方案差距不大. 可得在本研究區(qū)，物理相似法的效果隨著參證流域數(shù)量的增加而降低，與Viviroli等[48]發(fā)現(xiàn)的變化趨勢基本一致. 關(guān)注基于多流域物理相似法的方案3、6、9及方案4、7、10. 這些方案的模型效率損失均較大，且使用回歸法、物理估算法及其組合并不能有效提高多流域物理相似法的模型效率. 可見，使用其它方法代替物理相似法估算出部分參數(shù)，不會改變參證流域數(shù)量對參數(shù)估計效果的影響趨勢，不能有效彌補物理相似法參證流域數(shù)量選擇不當(dāng)造成的估計精度損失. 基于單流域物理相似法的方案效果最優(yōu)，下面聚焦于方案2、5、8. 在3個目標(biāo)流域，模型效率損失均是方案8<方案5<方案2，且在茅坪河與下牢溪，方案8的效率損失顯著低于方案2、5. 可見單獨使用物理相似法的參數(shù)估計效果最差；結(jié)合回歸法和物理相似法的參數(shù)估計效果與前者相比有小幅度提高；而根據(jù)模型參數(shù)特點綜合使用物理相似法、回歸法、物理估算法，可以顯著提高模型參數(shù)估計的精度. 說明在參證流域數(shù)量合適的情況下，本研究提出的綜合參數(shù)估計方案可以有效估計資料匱乏流域的新安江模型參數(shù).

考慮不同方案在研究區(qū)的整體效果及穩(wěn)定性，將3個目標(biāo)流域所有洪水場次的模型效率損失統(tǒng)計于圖4. 極少部分洪水場次存在估計參數(shù)略優(yōu)于率定參數(shù)的情況，因此這部分洪水的模型效率損失可能為負值. 由圖可見，方案2、5、8的模型效率損失中值最接近于0，且呈現(xiàn)明顯的負偏態(tài)分布，四分位距明顯小于其他方案. 說明相比于其他方案，方案2、5、8的參數(shù)估計效果更優(yōu)，且結(jié)果更穩(wěn)定. 而方案8的四分位距又顯著小于方案2、5，說明方案8的模型效率損失波動較小，估計參數(shù)的預(yù)報結(jié)果最為穩(wěn)定. 方案8在3個目標(biāo)流域的模型效率損失僅為0.26、0.28和0.01，基于洪水場次統(tǒng)計的模型效率損失中值僅為0.15，說明使用方案8估計的模型參數(shù)與率定的參數(shù)相比，精度僅有略微下降. 綜合評比后認為，方案8是最優(yōu)的模型參數(shù)估計方案，即單流域物理相似法結(jié)合回歸法及物理估算法是在本研究區(qū)內(nèi)進行新安江模型參數(shù)估計的最優(yōu)方法，可以用于該地區(qū)資料匱乏中小流域的洪水預(yù)報.

表5 不同參數(shù)估計方案結(jié)果統(tǒng)計

圖4 不同方案模型效率損失統(tǒng)計

5 結(jié)論及展望

將新安江模型應(yīng)用于長江中下游地區(qū)的32個山區(qū)中小流域，研究資料匱乏流域的參數(shù)估計方法. 本研究提出的參數(shù)估計方案包括物理相似法、回歸法、物理估算法及三者的不同組合，并同時考慮了參證流域的數(shù)量. 研究結(jié)果顯示：(1)在研究區(qū)內(nèi)，測量流域與目標(biāo)流域在空間上分為兩簇，空間分布具有顯著差異. 在區(qū)域化方法中，空間臨近法不適用，物理相似法適用，回歸方法只適用于估計參數(shù)CS及L. 物理估算法可以有效估計參數(shù)CS和KI；(2)目標(biāo)流域與測量流域空間距離較遠時，本研究提出的參數(shù)估計方案在不同程度上減少了參數(shù)估計導(dǎo)致的模型效率損失；(3)在本研究區(qū)，隨著參證流域數(shù)目的增加，物理相似法的估計效果下降. 使用回歸法、物理估算法及其組合并不能有效提升多流域物理相似法的模型效率，其它方法不能彌補物理相似法參證流域數(shù)量選擇不當(dāng)造成的模型效率損失；(4)單獨使用區(qū)域化方法往往不能得到理想效果. 在參證流域數(shù)量合適的前提下，物理估算法結(jié)合區(qū)域化方法，可以顯著提高模型參數(shù)估計效果；(5)提出的參數(shù)估計方案中，適用于長江中下游資料匱乏的山區(qū)中小流域的最佳方案為結(jié)合單流域物理相似法、回歸法、物理估算法的方案8.

測量流域與目標(biāo)流域空間距離較遠可以嘗試使用本研究提出的參數(shù)估計方案，但目前僅從主觀上認定其空間分布具有明顯差異，“較遠”的距離如何量化及空間距離是否存在適用性閾值等問題還需要進一步的研究. 同時，本研究中不同方案的參證流域數(shù)量僅選取了1、2、5個作為代表，多流域物理相似法的流域數(shù)目與參數(shù)估計效果的定量關(guān)系有待進一步研究. 本研究僅根據(jù)參數(shù)的物理意義估算出參數(shù)CS及L，新安江模型中的其它敏感參數(shù)能否通過類似方法進行估算，進一步提高資料匱乏地區(qū)的參數(shù)估計精度，還有待繼續(xù)研究.

6 附錄

附圖Ⅰ見電子版(DOI： 10.18307/2021.0223).