山洪預報模型歷史數據影響及解決方法探析

蘇曉林，秦根泉，黃少輝，林智偉，丁呈濤，3

（1.福州大學 土木工程學院，福建 福州 350000；2.四創科技有限公司，福建 福州 350000；3.福建省防災減災信息應用工程技術研究中心，福建 福州 350000）

0 引言

山區小流域洪水常在暴雨后發生，具有匯流時間短、水位陡漲陡落、成災時間快的特點，因而具有很強的突發性和破壞性，是我國防汛工作的難點和薄弱環節。隨著山洪災害防御體系在全國各地項目建設中陸續建立，我國山洪防治取得了許多成果。據統計，1961年以來，中國極端強降水事件呈增多趨勢，而2011年至2019年山洪災害年平均死亡人數近333人，較項目實施前的2000年到2010年的1179人大幅減少七成[1]。但是由于我國山區地形復雜、降雨時空分布不均，加上城鎮化的發展，導致山洪風險依舊嚴重。據統計，絕大部分已建的山洪監測預警系統還處于初級階段，隨著我國山洪災害防治非工程措施建設的深入，山洪監測預警系統建設正面臨從“感知型”向“智能型”的轉變，即從山洪災害調查評價得到的靜態風險區成果進行經驗性預報，向以山洪模型模擬成果進行預警轉變，以期得到更加精準和實時的山洪災害預測預警來指導避險[2]。

鑒于歷史水文資料的缺失錯漏是制約山洪預報模型成功搭建的重要因素，本文以海南白沙河流域為例，探索了一種通過選取典型流域進行洪水預報模型參數率定，再采用水文模型參數物理特性相似法將得到的參數方案移植到其它子流域的方法，以期解決山洪小流域洪水預報模型當前普遍存在的數據缺失問題。

1 研究現狀

1.1 山洪預報模型研究成果

當前我國的山洪預警預報方法總體上可分為兩種：一種是以統計歸納所得的臨界預警指標如臨界雨量閾值進行預報的經驗性方法；另一種是以模擬降雨產流過程的洪水預報模型進行預報的機理性方法[3]。隨著研究的發展，也出現了結合兩種方法的動態臨界預警指標法等。一般洪水預報模型從實現原理上可分為水文模型、水動力模型和人工智能預測模型。其中由于新興的人工智能預測模型的預測精度和輸入的歷史數據量有關，所以難以在普遍資料不足的山洪流域滿足精度要求。同樣的，嚴格基于物理過程的水動力模型具有更高的準確性，但因其對河道地形等資料的高要求，也不適合應用在無資料的山洪流域，而半經驗半物理過程的水文模型更加簡單、高效，在山區暴雨洪水模擬上應用更為廣泛[4]。

當前常用的洪水模擬水文模型方法主要有：新安江模型、HEC-HMS（Hydrologic Engineering Center-Hydrologic Modeling System） 模 型、TOPMODEL（Topgraphy Based Hydro-logical Model） 模 型、FFMS（Flash Flood Modul Simulation System）模型、NAM模型等[5-7]。新安江模型是趙人俊團隊研發的集總式水文模型，在較濕潤地區有較好的適用性，研究應用成果較多；HEC-HMS模型是美國陸軍工程兵團開發的半分布式次降雨徑流模型，在長時間尺度的模擬應用上較少，主要應用于濕潤和半濕潤地區；TOPMODEL模型是基于地形的半分布式流域水文模型，它考慮了地形因素，但忽略了土壤、植被等要素，主要適用于小流域尺度和半干旱地區；FFMS模型是劉昌軍團隊研發的模塊化小流域分布式水文模型，它是基于超滲、蓄滿時空動態組合的原理，主要運用于小流域暴雨洪水的模擬；DHI（Danish Hydraulic Institute）開發的集總式概念性的NAM模型，具有較簡單的結構和明確物理意義的參數，對國內外許多流域的降雨徑流過程都有較好的模擬效果，此外由于NAM模型對數據量的要求相對較少，因此在資料條件較差的區域也具有一定的可行性[8-11]。

1.2 模型數據存在的問題

依據相關規范在洪水預報中，必須使用至少10年的水文氣象資料，且還應包括各種代表性年份和一定場次洪水資料，以此作為精度評定，確定預報等級，其中對于驟發性洪水（如短歷時、高強度暴雨形成的小范圍山洪）的預報可采用氣象與水文相結合的方法，可不進行精度評定[12]。山洪預報雖然可以不進行精度評定，但在建立水文模型過程中，仍然需要較好的水文資料進行參數率定，才能得到反映地區情況較準確的水文模型參數方案。并且隨著監測手段和山洪治理建設的推進，對防洪安全和模型預報精度的要求也將有所提高。

雖然我國對山洪的治理工作已經取得初步成效，但由于山洪高發地區通常在較為偏遠的山區，限于較差的實際條件，通常難以得到足夠多的符合質量要求的歷史序列數據資料，使得山洪模型在構建和參數率定等方面存在許多困難。當前監測數據主要存在的問題有：①水文監測種類較少，河道通常只有水位數據缺乏相應斷面的水位流量關系；②水文資料的代表性不強，不能很好地代表流域情況；③水文資料存在不連續、不穩定的問題，數據質量不佳，降低了數據可靠性，也影響數據的篩選使用；④水文資料時間序列過少，數據量不足，難以滿足預報要求。

分析現有的數據問題，可以發現影響數據質量的原因主要有：①山區洪水流速快且含沙量大難以準確監測，監測設備不夠齊全，監測手段較為落后；②監測站點位置不夠合理，例如雨量監測站的密度不夠，不能準確代表整個流域的情況；③監測站點級別較低，多為自建站點，監測儀器的穩定性不足，存在缺測、誤測等問題；④山洪區域監測站點建設時間晚，監測時間短，缺乏多年長序列實測資料。

1.3 數據缺陷解決思路

在缺乏資料或資料質量欠佳的山洪高發地區，水文模型參數受限于現狀條件而難以率定，對此在實際工程中通常是采用水文模型參數區域化的方法，將其它有資料地區的水文模型參數參照運用于本地區中。當前水文模型參數區域化方法主要可以分為兩類：一類是參數移植法，主要有空間鄰近法和物理特性相似法等；另一類是參數回歸法，主要有多元回歸法和全局平均法等[13]。空間鄰近法是借用鄰近地區的模型參數的方法，它假定相近地區具有相似的水文條件，以空間距離作為主要判定依據，對于地形地貌變化不大的地區有較好的準確性。物理特性相似法是將地理條件、氣象條件等流域特征作為參數移用的判斷依據，與特征參數的選擇有較密切的關系。多元回歸法是建立模型參數與選取的地形特征參數的回歸方程，對特征參數選取和回歸方式有較高的要求，且需要較多的流域樣本。全局平均法是將流域內劃分的有資料的小流域所得參數進行算術平均，用以代表整個流域特征情況，一般適用于不大的流域。

由于參數回歸法對數據條件要求較高，因此當前應用更為廣泛的是參數移植法，而山洪預報機理性方法是基于水循環過程對降雨徑流過程的模擬，各基本假定和參數具有一定的物理意義，采用物理特性相似法更為合理。當前對流域相似性的判斷主要考慮氣候條件和下墊面特征[14]。對于處在同一流域中特別是同一小流域中的子流域來說，其氣候條件和下墊面特征通常都不會有太大的變化，具有較大的相似程度，在參數的移植借鑒上有著較高的可信度。NAM模型的參數在物理意義上主要與下墊面特征有關，對于具有相似下墊面特征的山洪區域，模型參數具有較小的差異性，在參數移植上有較好的可行性。

2 參數確定方法

2.1 研究區域概況

本文選取海南省的樂東縣白沙河流域開展研究。白沙河位于海南省的西南面，河道比降9.95‰。干流全長26km，流域集雨面積170km2。流域處亞熱帶海洋性季風氣候區，雨水十分充沛。白沙河流域中的水文站點數量較少，水位站點僅有兩個，一是流域出口處的丹村水位站（已損壞且缺乏歷史水位資料），另一個是支流上的牛尼頭水位站；而雨量站的數量也較少，且多分布于白沙河上支流區域。因此，對于整個白沙河流域洪水預報模型的搭建，站點分布及數據均不能滿足要求。子流域劃分采用ArcGIS水文分析工具處理白沙河流域DEM數據，以重要斷面（如水位站點）和支流匯入斷面作為流域劃分依據，將白沙河流域劃分為10個子流域。選取白沙河流域內歷史數據情況相對較好的牛尼頭水位站點控制的子流域5作為典型區域，子流域5內有尖峰嶺和牛尼頭兩個雨量站，根據子流域形狀和站點分布，由泰森多邊形作圖法確定其權重分別為0.9和0.1（圖1）。

圖1 白沙河流域劃分及站點分布圖Fig.1 Map of Baisha River Basin and hydrological stations

通過NAM模型對子流域5進行山洪模擬和參數率定，再通過對子流域典型特征的比對將參數移植到其它子流域模型。牛尼頭水位站建站年月為2013年1月，其上游控制子流域面積約為11km2。牛尼頭水位站的實測水位數據序列為2015年6月1日—2020年12月8日，有5.5年，本次基于水位數據變化和降雨數據變化相匹配的原則，對數據準確性和完整性進行判斷，篩選可用洪水場次，結果詳見表1。

表1 牛尼頭水位站水位監測數據情況表

從表1可知，該站點的實測水位序列不足，質量也欠佳，中間存在許多缺測和誤測數據，而且雨量站存在部分降雨場次的缺失，綜合考慮歷史資料的完整性和代表性后選取3場較為可信的降雨和水位數據供模型參數率定和驗證。

2.2 參數率定方法

2.2.1 模型原理

NAM模型將降雨徑流過程分為四個相互影響的儲水層的計算，其具體結構見圖2[15]。

圖2 NAM模型計算結構圖Fig.2 NAM model calculation structure diagram

模型的參數是對流域平均情況的代表取值，一般無法通過實測獲得，只能通過率定得到適合該流域的參數方案，參數物理意義及取值范圍見表2[16]。

表2 NAM模型主要參數意義及取值范圍

2.2.2 模型率定過程

水文產匯流模型參數率定需要將實測水位數據轉換成對應流量數據。根據實測數據，牛尼頭監測斷面屬于寬淺型河道，且斷面附近河道較為順直，斷面尺寸變化不大，可以采用如式（1）的曼寧公式計算水位流量關系。

式中，Q為流量，單位m3·s-1；A為斷面過流面積，單位m2；v為過流斷面平均流速，單位m·s-1；R為水力半徑；n為糙率，取值0.040；J為洪水水面線比降，由于缺乏資料，實際采用河道比降，取值8.8‰。

模型參數的確定，采用自動率定和人工率定相結合的方式。

自動率定采用SCE-UA隨機搜索參數優化算法，該算法自1992年由Duan提出以來在水文模型參數率定中得到廣泛應用和驗證[17]。自動率定的目標為：①總水量平衡，即模擬的徑流量與實測值相等；②洪水過程線的形狀相近，變化規律類似；③洪水流量的峰值接近，主要考慮洪峰出現時間、洪峰流量大小以及水量等條件；④低流量吻合，即基流情況吻合[15]。

在進行自動率定之后，再進行人工率定，反復進行參數調整步驟，最終得到最符合各場次洪水的總體情況的參數方案。人工調整步驟如下：①通過調整參數Lmax和Umax來平衡水量；②通過調整參數CK1、CK2和CQOF來調整洪峰流量的大小；③通過調整參數CKBF來調整基流[15]。使用Nash-Sutcliffe系數評價模擬結果與實測資料的符合程度，計算公式見公式2。

式中，R2為Nash-Sutcliffe系數，越接近1表示模擬結果和實測資料的差異性越小；Qobs（i）為實測流量，單位m3·s-1；Qsim（i）為模擬流量，單位m3·s-1；為實測流量平均值，單位m3·s-1；N為實測流量資料序列的長度；wi為權重系數。

3 參數結果與應用

3.1 參數率定結果分析

根據站點現有數據情況，盡量多選取可用于模型計算的山洪場次進行模型率定和驗證。經分析，選取2016年的2場山洪進行模型參數率定，2017年的1場山洪作為驗證，參數率定的模型模擬結果見表3和圖3-4。基于同一模型參數方案，不同場次降雨間模擬效果差異的原因還可為山洪預報對降雨分布的敏感性，不同場次降雨雨強中心和站點斷面的距離存在差異，這就導致雨量站的實際分配效果和由泰森多邊形計算出的雨量分配權重方案有所差異，且將直接影響到匯流時間的長短，從而引起模型參數之間的差異。同時降雨強度的不同會影響徑流量和下滲量的分配比例，也是影響模型參數的原因之一。

從表3中可以看出，洪水率定期選取的2個場次洪水的模擬洪峰相對誤差均在洪峰許可誤差20%之內，符合洪峰模擬精度要求。從圖3、4的模擬結果中可以看出，模擬洪水過程和實測情況整體吻合性較好，基本上能夠正確地反映出洪水過程，出現洪峰的時間也較為接近。一般情況下洪峰出現時間要晚于降雨峰值出現的時間，這一雨洪滯后性也在模擬結果中得到了驗證。

表3 牛尼頭站斷面NAM模型場次山洪模擬結果

圖3 20160924場次洪水實測和模擬結果對比圖Fig.3 Comparison of measured and simulated results of 20160924 floods

圖4 20161018場次洪水實測和模擬結果對比圖Fig.4 Comparison of measured and simulated results of 20161018 floods

但是，其中20160924場次洪水模擬的Nash-Sutcliffe系數不佳，洪水模擬過程也不甚理想。分析原因是在率定參數的過程中，為了使2場洪水總體模擬效果較好，在自動率定之后進行手動參數調整。考慮到20161018場次洪峰流量更大，風險程度更大，所以在參數調整過程中，更加側重該場次的模擬效果，在保證洪峰誤差和峰現時間誤差滿足要求的前提下，適度犧牲20160924場次模擬精度。

受限于較少的可用洪水場次，且考慮到偶然性的影響，難以給出多種工況情景下參數方案，同時也為了便于模型的部署應用，綜合考慮現有數據資料率定結果，最終確定一套參數方案，如表4所示，各參數均在合理范圍之內，將被作為該地區的參考方案。

表4 山洪模擬NAM模型主要參數取值方案

3.2 模型模擬驗證分析

選取2017年的一場山洪進行參數驗證，模擬結果見表5和圖5。從表5可以看出，驗證期的洪水場次模擬洪峰誤差和峰現時間誤差基本滿足要求，Nash-Sutcliffe系數滿足要求，此外從圖5的模擬結果也可以看出洪水過程的模擬情況同實測較為吻合，說明采用的參數方案具備一定的準確性。

圖5 20170826場次洪水實測和模擬結果對比圖Fig.5 Comparison of measured and simulated results of 20170826 floods

表5 牛尼頭站斷面NAM模型場次山洪模擬驗證結果

3.3 參數移植方法

3.3.1 參數移植應用思路

由于不同山洪流域的土壤和地形可能存在較大差異，因而不能直接使用代表性區域的參數方案，本次研究根據物理特性相似法制定NAM模型在山洪流域的參數移植方法，參數移植思路分兩種：一是分析模型參數的敏感性，篩選對模擬結果影響較大的模型參數和相關的流域物理特征，并建立模型參數和對應流域物理特征之間的數學關系，以流域物理特征計算模型參數取值；二是選取足夠數量的具有不同流域物理特征的研究區分別進行參數率定，進而組合成一套適用多種類型流域的參數方案，建模時通過流域物理特征選取參數方案。本研究區域流域物理特征如表6所示。

表6 白沙河牛尼頭小流域物理特征表

要建立整個白沙河流域的降雨徑流模型，需要先建立各個子流域的水文產匯流模型，然后建立各個子流域匯流出口到流域總出口斷面的河道水動力模型。由于在白沙河流域內，除牛尼頭子流域外的其它子流域都缺乏站點和資料，采用上述方法二，以牛尼頭子流域的參數方案為基礎，對模型參數進行調整，移植運用到其它子流域，從而搭建出整個白沙河流域的降雨徑流模型。

3.3.2 參數敏感性與相關流域物理特征分析

模型參數的敏感性指的是參數變化對模擬結果影響的大小，在現有的NAM模型的參數敏感性分析研究中，姚璇[18]采用基于擾動分析法的相對靈敏度方法對NAM模型9個主要參數進行了局部敏感度分析，結果表明參數Lmax、Umax、CQOF和TG（參數含義參見表2）具有較高的敏感性，而CK1和CK2在日徑流模擬結果影響較小。陳福容等[19]則認為參數Lmax、CK1和CK2是較為敏感的參數。此外李磊等[20]研究發現參數Lmax、Umax、CQOF對徑流總量影響明顯，CK1、CK2對峰值的影響明顯。研究普遍認為參數Lmax、Umax、CQOF是對水量平衡有較大影響的參數，而影響匯流時間的CK1、CK2同樣具有較大的敏感性，且在本次手動率定過程中同樣發現了類似的規律，因此對模擬結果影響較大的這5個參數，需要在參數移植中重點考慮。根據各參數的物理意義，總結其主要相關的流域特征參數，如表7所示。

表7 NAM模型部分參數及主要相關流域特征參數表

參數Lmax、Umax分別是淺層蓄水層最大含水量和地表蓄水層最大含水量，通常前者約為后者的10倍，主要與土壤有關，且強烈受到植被的影響。與洪峰流量有顯著相關的CK1、CK2和CQOF這三個參數和流域特征之間一般具有如下關系：當土壤透水率越大時，CK1、CK2和CQOF越小；當坡度越大時，CK1、CK2越小，而CQOF越大；當植被覆蓋率越大時，CK1、CK2也越大，而CQOF越小；當集水面積越大而流域形狀系數越小時，CK1、CK2越大，而CQOF變化較小。

3.3.3 白沙河子流域參數調整取值

白沙河各子流域的流域特征情況基本類似，因此，可以考慮在牛尼頭子流域模型參數方案的基礎上調整Lmax、Umax、CQOF、CK1和CK2參數，而保持其它參數不變。白沙河各個子流域的模型參數取值方案如表8所示。

表8 白沙河各子流域物理特征和模型參數取值表

采用上述參數方案建立白沙河各子流域的NAM模型，得到各子流域匯流的徑流量，再通過水動力模型將其匯流到流域總出口斷面（即子流域9匯流斷面處），得到流域總出口斷面的徑流量過程。圖6為干流沿程各斷面位置處的流量模擬結果。從圖6可以看出，從河源到流域出口，沿程河道斷面的洪水流量不斷增加，且峰現時間也逐漸延后，符合洪水演進的正常規律，后期將結合流域出口新建的水情站點監測數據，對模型參數進行優化，逐步提高模型預報精度。

圖6 白沙河干流沿程流量模擬結果圖Fig.6 Flow simulation results of Baisha River mainstream

4 結論與展望

本文主要探討了山洪小流域地區的歷史資料數量不足和質量欠佳情況下洪水預報模型的構建過程。針對山洪流域下墊面變化劇烈，建設站點數量不足的現狀，本文提出區域化代表參數的思路，并采用NAM模型進行模擬研究，是一種基于現狀，解決實際問題的建模方法探索。

（1）為解決山洪預報模型歷史數據質量的限制，本文以白沙河流域為試點，選取數據條件相對較好的牛尼頭站點控制的子流域為代表進行山洪預報模型參數的率定和驗證，根據山洪流域參數移植方法的適用性選取物理特征相似法進行參數移植，從中得出一套較為合理的參數方案供該區域的山洪預報模型建模使用。

（2）本研究運用基本的技術手段，進行模型參數的率定驗證和移植，以牛尼頭子流域率定參數方案作為其它子流域參數方案的調整基礎，可以避免其它子流域參數取值的盲目性，在無資料驗證的情況下，最大限度地提高模擬的準確性。

（3）對于耦合水文模型和水動力模型的降雨徑流模型，總出口控制斷面的模擬精度是氣象預報、水文模擬和水動力模擬三個階段的精確性的最終體現，在模型模擬精度分析時，需要對各個過程綜合考慮來提高模擬預測的精度。