李 姣,王麗榮,王 潔,鄧 鵬,趙鐵松

(1. 河北省氣象災(zāi)害防御和環(huán)境氣象中心,河北 石家莊 050021; 2. 河北省氣象與生態(tài)環(huán)境重點實驗室,河北 石家莊 050021;3. 南京信息工程大學(xué) 水文與水資源工程學(xué)院,江蘇 南京 210044)

0 引言

山洪災(zāi)害具有來勢猛、預(yù)見期短、危害性大、預(yù)報困難等特點。我國山丘區(qū)面積大、承災(zāi)體脆弱,一旦發(fā)生山洪災(zāi)害,往往會造成嚴(yán)重的經(jīng)濟損失以及人員傷亡[1],山洪預(yù)警預(yù)報相關(guān)研究的開展,對降低山洪災(zāi)害危險性,具有現(xiàn)實應(yīng)用價值。

現(xiàn)行的山洪災(zāi)害預(yù)警主要依據(jù)臨界雨量,即當(dāng)降雨量超過某一臨界值就會引發(fā)山洪災(zāi)害,如何準(zhǔn)確地估算臨界雨量在山洪災(zāi)害預(yù)警預(yù)報工作中至關(guān)重要[2]。國內(nèi)外學(xué)者在臨界雨量確定方法方面提出了不同思路[3-5],如美國國家氣象局水文研究中心,基于典型的降雨、產(chǎn)流匯流等,考慮了降雨、流域下墊面特性等因素,利用Flash Flood Guidance(FFG)[3]反推臨界雨量。國內(nèi)山洪臨界雨量的確定主要采用統(tǒng)計方法、物理計算和水文模擬法等,陳桂亞等[6]、樊建勇等[7]、王仁喬等[8]分別以湖南、江西、湖北小流域為例,對山洪災(zāi)害過程的降雨數(shù)據(jù)進行統(tǒng)計,分析得到了研究區(qū)山洪災(zāi)害臨界雨量;江錦紅等[9]計算山洪災(zāi)害前期累計降雨量和災(zāi)前1 h的降雨量,繪制暴雨臨界曲線,判定是否達到山洪災(zāi)害臨界雨量;李磊[10]、王新宏等[11]基于暴雨與洪水同頻率假定,采用水位/流量反推法計算無資料區(qū)/資料匱乏地區(qū)的山洪災(zāi)害臨界雨量;劉淑雅等[12]利用新安江模型模擬了撫河流域不同時段、不同初始土壤條件下的臨界雨量;羅倩等[13]采用HEC-HMS水文模型試算湖北省高堤河流域的臨界雨量,并分析了不同的降雨空間分布對臨界雨量計算的影響。以上國內(nèi)臨界雨量計算方法得到的多為靜態(tài)雨量,而實際情況下,受下墊面狀況、前期降雨、土壤濕度、河道流量等的影響,臨界雨量是動態(tài)變化的。近年來,水文氣象學(xué)家基于水文模型[14-16]研究開展了動態(tài)臨界雨量的研究,如劉志雨等[14]以GBHM水文模型為基礎(chǔ),計算了不同土壤飽和度下的不同歷時的江西遂川江流域山洪災(zāi)害動態(tài)臨界雨量,但是相關(guān)研究較少,且多數(shù)研究僅考慮了前期土壤含水量和降雨量的影響,未考慮前期河道流量和產(chǎn)流面積因素,且南方地區(qū)研究較多,北方山丘區(qū)關(guān)于山洪動態(tài)臨界雨量的研究較少。

本文選取位于河北太行山區(qū)的坡底小流域,基于新安江模型,綜合考慮了流域前期的土壤飽和度、前期河道流量和產(chǎn)流面積,構(gòu)建了山洪災(zāi)害動態(tài)臨界雨量計算模型,可逐小時計算未來1、3、6、12、24 h這5種不同歷時的I-IV級的動態(tài)臨界雨量,結(jié)合氣象部門逐時降水預(yù)報,可實現(xiàn)對山洪災(zāi)害的逐時滾動預(yù)報,對山洪業(yè)務(wù)提供技術(shù)支撐。

1 研究流域與數(shù)據(jù)資料

1.1 研究區(qū)概況

選取河北省太行山區(qū)坡底流域為試點,進行山洪災(zāi)害動態(tài)臨界雨量閾值試算。坡底流域隸屬于海河流域,流域面積283 km2;年降水量在400～700 mm之間,屬于半濕潤地區(qū),降水量年內(nèi)變化大,主要集中在6—9月;流域地形以山地為主,海拔305～1742 m,地勢西高東低、坡度變化較大;流域植被覆蓋較好;坡底流域共有9個降水觀測站和1個水文站(圖1)。受其陡峭的地形及氣候特征影響,該流域的山洪、洪水災(zāi)害頻發(fā)[17],出現(xiàn)過“63·8”、“96·8”、“16·7”等特大洪水災(zāi)害。

圖1 坡底流域海拔和站點分布

1.2 研究資料

坡底流域1996、2000、2016、2018年汛期(7—8月)洪水過程的逐小時流量數(shù)據(jù),由邢臺水文勘測研究中心提供;與山洪過程時段對應(yīng)的流域內(nèi)9個雨量站的逐小時降水?dāng)?shù)據(jù),以及離坡底流域最近的邢臺氣象站的逐小時蒸發(fā)數(shù)據(jù),降水和蒸發(fā)數(shù)據(jù)由河北省氣象局提供;所有數(shù)據(jù)已經(jīng)經(jīng)過質(zhì)量控制,數(shù)據(jù)可靠,分布均勻(圖1)。

2 研究方法

2.1 新安江模型

新安江模型已被廣泛應(yīng)用于濕潤與半濕潤地區(qū)的山洪過程模擬[18-19],模型結(jié)構(gòu)分為蒸散發(fā)、產(chǎn)流、分水源和匯流4個層次進行計算:采用3層模型計算流域蒸散發(fā)量;按蓄滿產(chǎn)流概念計算流域產(chǎn)流量;將水源劃分為飽和地面徑流、壤中水徑流和地下水徑流三類;針對3種徑流類型,分別選取單位線法(地面徑流)和線性水庫法(壤中水和地下水徑流)進行匯流計算,河道匯流采用分段連續(xù)演算的馬斯京根算法[20-21]。

模型的輸入數(shù)據(jù)包括流域面雨量、蒸發(fā)量、流域特征參數(shù)、模型參數(shù)等,輸出數(shù)據(jù)為流域的土壤飽和度、河道流量、產(chǎn)流面積等,利用山洪過程率定新安江模型參數(shù)。

2.2 動態(tài)臨界雨量計算流程

本研究采用由臨界流量反推臨界雨量的方法。首先確定研究區(qū)域的山洪災(zāi)害臨界流量,根據(jù)確定的致災(zāi)臨界流量和當(dāng)前時刻的流域土壤飽和度、河道流量、產(chǎn)流面積等,查找可產(chǎn)生臨界流量的降水,作為臨界雨量。受到前期降水、蒸發(fā)等影響,不同時刻流域的土壤飽和度、河道流量和產(chǎn)流面積是動態(tài)變化的,因此,不同時次的山洪災(zāi)害臨界雨量也是動態(tài)變化的,即:動態(tài)臨界雨量。

1)臨界流量確定

根據(jù)流域出口的流量大小判定流域是否發(fā)生山洪。由于缺乏長時間序列的流量資料,本研究采用百分位數(shù)法確定不同等級的山洪災(zāi)害臨界流量值。將流域內(nèi)的流量數(shù)據(jù)按大小排序,取其85%、80%、75%、70%位數(shù)對應(yīng)的流量作為I-IV四個等級的山洪災(zāi)害臨界流量值,代表發(fā)生山洪災(zāi)害的可能性逐級減小。

H8：品牌影響力對新疆農(nóng)產(chǎn)品品牌競爭力有正向影響，即新疆農(nóng)產(chǎn)品品牌知名度、美譽度、忠誠度越高，新疆農(nóng)產(chǎn)品區(qū)域品牌競爭力越強。

2)動態(tài)臨界雨量確定

將計算時刻前≥15 d的降水、蒸發(fā)、臨界流量以及率定好的模型參數(shù)輸入新安江模型,模擬得到當(dāng)前時刻的土壤飽和度、河道流量和產(chǎn)流面積,在此基礎(chǔ)上,再采用搜索效率較高的二分法,搜尋能夠剛好達到山洪臨界流量對應(yīng)的降雨量,即為當(dāng)前時刻山洪災(zāi)害臨界雨量值(圖2),具體計算步驟如下:

圖2 動態(tài)臨界雨量模型計算流程

步驟1:給定降雨下限P1(例如:0.1 mm)和降雨上限P2(例如:1000.0 mm),求出其中位數(shù)P;

步驟2:將前期觀測降水?dāng)?shù)據(jù)和P作為新安江模型的降水輸入,模擬降雨過程中的最大流量Qm;

步驟3:將最大流量Qm與山洪災(zāi)害臨界流量Qc進行比較,如果Qm>Qc,將P替換掉P2,作為降水上限,反之,將P換掉P1,作為降水下限;

步驟4:重復(fù)步驟1-3,直到P2與P1的差值小于降水?dāng)?shù)據(jù)精度(一般精確到0.1 mm),那么,根據(jù)此時的雨量P模擬計算得出的流量最大值最接近山洪臨界流量,即,雨量P為Qc相對應(yīng)的臨界雨量。

3)臨界雨量時段及等級確定

為了滿足山洪災(zāi)害預(yù)警的需求,將預(yù)警時段確定為1、3、6、12、24 h這5種歷時。根據(jù)上述確定的山洪災(zāi)害臨界流量值,動態(tài)臨界雨量也分為I—IV 4個等級,代表發(fā)生山洪災(zāi)害的風(fēng)險依次降低。

3 結(jié)果與分析

3.1 水文模型參數(shù)率定

利用新安江模型開展山洪過程模擬,首先要進行模型參數(shù)的率定,利用1996、2000、2016、2018年間的11次山洪過程進行參數(shù)率定,并驗證模型參數(shù)。新安江模型參數(shù)物理意義明確,取值受地區(qū)影響較大。其中,有的參數(shù)很敏感,稍微變化就對輸出有很大影響;有的參數(shù)反應(yīng)遲鈍,對輸出的影響很小,參數(shù)敏感程度受濕潤程度、水量大小等影響。根據(jù)模型結(jié)構(gòu),按照蒸散發(fā)計算、產(chǎn)流與水源劃分、匯流演算的次序?qū)δＰ蛥?shù)進行分塊率定。根據(jù)流域面積、自然地理特征等,給出模型參數(shù)的初始值,運行水文模型,模擬逐小時河流斷面流量;對比實測與模擬的水量平衡、洪水過程線、洪峰流量、峰現(xiàn)時間、確定性系數(shù),不斷調(diào)整模型參數(shù),直至滿足預(yù)報精度要求[22]。表1為坡底流域模型參數(shù)率定結(jié)果。

表1 坡底流域模型參數(shù)率定結(jié)果

3.1.2 山洪過程模擬結(jié)果

表2給出了坡底流域11次山洪過程的模擬結(jié)果,包括率定期山洪7場,驗證期4場?？梢钥闯?對于平均徑流量模擬結(jié)果而言,有10場山洪的模擬相對誤差小于15%,僅有一場“16·7”大洪水接近20%;從洪峰流量的模擬結(jié)果來看,所有場次洪水的模擬結(jié)果較好,相對誤差均小于13.6%,峰現(xiàn)時間均控制在2 h以內(nèi);所有場次洪水的實測流量與模擬流量的確定性系數(shù)均高于0.67。根據(jù)水情預(yù)報規(guī)范[23],以確定性系數(shù)為評價標(biāo)準(zhǔn),確定性系數(shù)在0.5以上的合格預(yù)報次數(shù)為11次,合格率達100%,確定性系數(shù)在0.7以上(模擬準(zhǔn)確度達到乙級)的場次為8次,占總體的72.7%。

以“16·7”的洪水過程為例分析場次洪水模擬情況。2016年7月18日22時—20日10時,坡底流域出現(xiàn)持續(xù)性強降水,流域36 h累計面雨量達230.7 mm,圖3給出了此次山洪過程的降水實況以及過程流量和土壤飽和度的模擬結(jié)果。圖3(a)可以看出,模擬與實測流量過程基本一致,洪峰出現(xiàn)時間基本一致,但是,模擬洪峰流量(1520.8 m3/s)小于實測洪峰流量(1630 m3/s),這有可能是模型自身模擬誤差所致,也有可能是在使用雨量站數(shù)據(jù)作為模型輸入數(shù)據(jù)時,由于雨量站布設(shè)密度不夠,不能全面監(jiān)測全流域的空間降雨情況,導(dǎo)致模型模擬輸入雨量偏小,因此模擬流量偏小[24];圖3(b)可以看出,土壤在開始降水一段時間后達到飽和狀態(tài),降水停止2 d后,土壤飽和度緩慢降低。整體來看,本次洪水過程模擬基本符合實際情況,平均徑流量模擬相對誤差為19.8%,洪峰流量相對誤差為6.7%,峰現(xiàn)時差為1 h,確定性系數(shù)為0.67,符合洪水預(yù)報精度規(guī)定。

圖3 坡底小流域“16·7”洪水過程模擬結(jié)果

綜合以上分析,新安江模型在坡底流域的模擬結(jié)果較好,可以用于山洪動態(tài)臨界雨量的計算。

3.2 動態(tài)臨界雨量的估算

采用百分位數(shù)法,計算坡底流域流量數(shù)據(jù)的85%、80%、75%、70%位數(shù)的流量值,分別為357、139、71、58 m3/s,作為I—IV等級山洪災(zāi)害臨界流量。根據(jù)山洪代表性的原則,選取2016年7月19日特大山洪災(zāi)害過程前后的3個不同預(yù)報時刻以及2018年中小規(guī)模山洪的場次的預(yù)報時刻,計算了18個不同事件的山洪臨界雨量閾值。首先,將每個時刻的臨界雨量閾值與實際面雨量進行比對,來判定是否發(fā)布山洪災(zāi)害預(yù)警及預(yù)警等級;然后,將對應(yīng)時刻的臨界流量閾值與實際流量進行比對,判定實際是否發(fā)生山洪。

表3給出了計算結(jié)果,包括是否發(fā)生山洪災(zāi)害預(yù)警、山洪災(zāi)害預(yù)警等級等。從預(yù)報是否發(fā)生山洪來看,18次事件中,有15次預(yù)報準(zhǔn)確,1次空報,2次漏報,預(yù)報準(zhǔn)確率為83.3%;其中,漏報的事件分別屬于6、24 h長歷時的預(yù)警預(yù)報,雖然長歷時沒有預(yù)報準(zhǔn)確,但其前面的短歷時預(yù)報卻準(zhǔn)確及時地對山洪發(fā)出了預(yù)警預(yù)報。從山洪預(yù)警等級來看,15次成功預(yù)報山洪災(zāi)害事件中,有11次準(zhǔn)確預(yù)報,4次預(yù)報等級偏高,偏高時次為7月19日01時的12 h預(yù)報、19日14時的1 h預(yù)報和8月24日20時的3、6 h預(yù)報,另外,8月24日20時1 h出現(xiàn)空報,可以看出,預(yù)報偏高時刻主要出現(xiàn)在降雨過程開始或者雨強開始增強的階段,這可能與水文模型的模擬誤差有關(guān),模擬流量大于實際流量,導(dǎo)致預(yù)報等級偏高,而水文模型模擬效果不理想的原因很多,其中最主要的原因是山洪場次數(shù)據(jù)偏少,樣本代表性較差,導(dǎo)致率定的模型參數(shù)存在誤差。

表3 山洪災(zāi)害預(yù)報效果統(tǒng)計

圖4給出了“16·7”暴雨過程前后坡底流域不同時刻的山洪災(zāi)害臨界雨量,可以看出,受到前期降水的影響,山洪災(zāi)害臨界雨量是隨時間動態(tài)變化的。以表3中選取時次的1 h時段的臨界雨量計算結(jié)果為例進行分析,7月19日01時,山洪災(zāi)害發(fā)生前,1 h的I—IV級臨界雨量閾值分別為30、16、11、10 mm;隨著降水的增大,山洪臨界雨量閾值逐漸減小,至19日14時,1 h的臨界雨量閾值變?yōu)?、0、0、0 mm,代表該地區(qū)正在遭受山洪災(zāi)害;山洪過程結(jié)束之后,不同歷時的山洪災(zāi)害臨界雨量閾值逐漸增大,8月24日再次出現(xiàn)降水,至24日20時,1 h的臨界雨量閾值分別為27、11、3、2 mm(表3),當(dāng)前時刻的1 h臨界雨量閾值介于7月19日01時與14時之間,這是由于受到前期降水影響,流域土壤飽和度、河道流量、產(chǎn)流面積發(fā)生改變的結(jié)果。

圖4 坡底流域“16·7”洪水過程動態(tài)臨界雨量閾值時間序列圖

以7月19日01時的臨界雨量閾值模擬結(jié)果看某一時刻的閾值特征,可以看出,臨界雨量大小是隨著不同時間尺度變化的,隨著歷時與預(yù)警等級的增高,臨界雨量逐漸變大,這種變化并不是線性的(圖5),這說明導(dǎo)致山洪災(zāi)害發(fā)生所需要的降雨量并不是與降雨歷時成倍數(shù)增加的,歷時越短,導(dǎo)致山洪發(fā)生所需的雨強越大,即1 h發(fā)生山洪的雨強最大,24 h發(fā)生山洪的雨強最小,這與流域本身的調(diào)蓄作用有關(guān)。

圖5 坡底流域2016年7月19日01時山洪災(zāi)害臨界雨量閾值

通過以上分析可見,模型整體預(yù)報準(zhǔn)確率尚可,計算得到的山洪災(zāi)害臨界雨量閾值隨時間和降水量的變化趨勢正確,相比靜態(tài)臨界雨量,更能準(zhǔn)確地反映出計算時刻流域面臨的山洪災(zāi)害風(fēng)險大小。

4 結(jié)論與討論

以往的山洪業(yè)務(wù)中用到的致災(zāi)臨界雨量閾值多為靜態(tài)閾值,但是受到前期降雨影響,土壤飽和度、河道流量和產(chǎn)流面積都是動態(tài)變化的,靜態(tài)閾值將不滿足于業(yè)務(wù)需求,本文以新安江模型為基礎(chǔ),開發(fā)了動態(tài)山洪臨界雨量計算模型,綜合考慮前期影響因素,可計算多個時間尺度多等級的山洪動態(tài)臨界雨量閾值,不僅為坡底流域山洪預(yù)警業(yè)務(wù)提供技術(shù)保障,也為其他流域山洪預(yù)警提供了研究思路,主要研究結(jié)論如下:

1)新安江模型可以很好地模擬出坡底流域的洪水過程,確定性系數(shù)在0.67以上,峰現(xiàn)時間控制在2 h以內(nèi),平均徑流量相對誤差和洪峰流量模擬誤差均在合理范圍內(nèi)。

2)計算得到的動態(tài)臨界雨量隨降雨時段不同動態(tài)變化,其變化強度受雨強、累計雨量、降雨時長等影響?？傮w來講,雨強越強、累計雨量越大,臨界雨量閾值變化越明顯。

3)導(dǎo)致山洪災(zāi)害發(fā)生所需要的降雨量并不是與降雨歷時成倍數(shù)增加的,歷時越短,導(dǎo)致山洪發(fā)生所需的雨強越大,即1 h發(fā)生山洪的雨強最大,24 h發(fā)生山洪的雨強最小。

4)本研究構(gòu)建的山洪災(zāi)害動態(tài)臨界雨量計算模型,對判別是否發(fā)生山洪災(zāi)害,預(yù)報準(zhǔn)確率可達83.3%,預(yù)警效果較好。