劉昌軍,文 磊,周 劍,趙懸濤,郭 良,魏永強
(1.中國水利水電科學(xué)研究院 減災(zāi)中心,北京 100038;2.河海大學(xué) 水文水資源學(xué)院,江蘇 南京 210098;3.中國科學(xué)院 西北生態(tài)環(huán)境資源研究院,甘肅 蘭州 730000;4.華北水利水電大學(xué) 水利學(xué)院,河南 鄭州 450045;5.湖南省水利水電科學(xué)研究院,湖南 長沙 410007)
我國山洪災(zāi)害點多面廣、突發(fā)性強、破壞力大,且多發(fā)生在偏遠山區(qū)、交通不便、通訊不暢。山洪災(zāi)害是嚴重威脅群眾生命財產(chǎn)安全的主要災(zāi)種之一,日漸成為防洪減災(zāi)的短板,防御工作一直是我國防汛和防災(zāi)減災(zāi)工作中的難點和薄弱環(huán)節(jié)[1]。近年來無資料小流域暴雨山洪模擬研究得到快速發(fā)展[2],目前常用的是采用水文學(xué)方法和水動力學(xué)方法[3]?;诜植际剿哪P偷男×饔虮┯旰樗M得到快速發(fā)展和廣泛應(yīng)用[4-5]。黃燕等[6]開展了水文模型在山洪模擬中的比較應(yīng)用研究,研究結(jié)果發(fā)現(xiàn)現(xiàn)有水文模型在模擬短歷時山洪效果一般。
隨著計算機技術(shù)的快速發(fā)展,采用水動力學(xué)方法進行暴雨山洪過程的模擬是今后研究方向和熱點[7-8]。平面二維數(shù)值模型已成為洪水演進模擬的主流方法,但間斷水流、干濕邊界處理和計算效率仍存在較多難點,是急需解決的關(guān)鍵問題[9]。潘佳佳等[10]認為基于運動波和擴散波等簡化水動力學(xué)模型因忽略了流動的力學(xué)因子,不能較好的模擬山洪的形成與演化機理。張培文等[11]采用有限元法研究了降雨條件下坡面徑流和入滲問題的耦合模擬;王嬌進行了坡面徑流—入滲耦合的解析解及無網(wǎng)格數(shù)值模擬研究[12];基于Godunov格式的以求解黎曼近似解的淺水數(shù)值模擬技術(shù)已逐漸成熟,該方法能模擬光滑解,又能自動適應(yīng)復(fù)雜流態(tài)變化,并得到廣泛應(yīng)用[13-14]。張大偉等[15]提出了新的基于Godunov格式的地表徑流二維水動力模型,使干濕邊界問題得到很好解決。張奕虹等[16]提出了新的水深重構(gòu)方法,解決了復(fù)雜地形上暴雨山洪數(shù)值模擬常面臨的計算水深小于零、非物理大流速等問題。闞光遠等[17-18]在GPU并行算法方面作了大量研究,發(fā)現(xiàn)使用GPU并行算法可大幅提高計算效率。GPU、CPU等并行計算技術(shù)已經(jīng)在水文和水動力學(xué)技術(shù)中得到初步應(yīng)用,具有廣闊的應(yīng)用前景[19]。
針對暴雨山洪模擬計算問題,以湖南瀏陽市寶蓋寺小流域2012年降雨徑流過程為例,本文采用筆者提出的時空變源混合產(chǎn)流模型和模塊化小流域暴雨山洪分布式水文模擬模型以及基于Godunov格式的二維坡面徑流數(shù)值模擬模型對小流域暴雨山洪進行了對比研究,從建模過程、計算效率和應(yīng)用等方面,對比分析了兩種方法計算精度和計算效率及計算結(jié)果合理性,系統(tǒng)分析了坡面產(chǎn)匯流機制和微地形對降雨徑流過程影響,為無資料小流域暴雨洪水計算和預(yù)報預(yù)警提供技術(shù)借鑒。
在山洪災(zāi)害調(diào)查評價成果和實時水文氣象大數(shù)據(jù)基礎(chǔ)上,筆者提出了以實時衛(wèi)星、雷達和臺站等多源降水融合數(shù)據(jù)和預(yù)報降雨數(shù)據(jù)為驅(qū)動,以調(diào)查評價小流域拓撲關(guān)系及基礎(chǔ)屬性為基礎(chǔ),以7類水文單元(流域、河段、節(jié)點、分水、水源、洼地、水庫)水力關(guān)系為核心,以歷史暴雨洪水和率定參數(shù)為先驗知識,以機器學(xué)習(xí)和并行計算為手段,以降雨蒸發(fā)、產(chǎn)匯流、河道演進、水庫調(diào)蓄等水文過程為主線,基于人工智能和大數(shù)據(jù)技術(shù),構(gòu)建以小流域為計算單元,集成模型庫、人工智能算法庫、參數(shù)庫和先驗知識庫為一體的模塊化分布式水文模型(FFMS)[20],模型計算流程和結(jié)構(gòu)見圖1。
圖1 分布式水文模型流程和結(jié)構(gòu)
2.1 時空變源混合產(chǎn)流模型圍繞山丘區(qū)中小流域地形地貌多樣、產(chǎn)流機制時空變化復(fù)雜等問題,劃分了山丘區(qū)山坡地貌水文響應(yīng)單元,采用基于一維入滲理論的包氣帶土壤非線性下滲計算方法,計算濕潤鋒下移過程及入滲量,提出了超滲、蓄滿時段轉(zhuǎn)換準則,建立了在平面、垂向、時段上時空變源混合產(chǎn)流模型[21],以精確模擬山丘區(qū)中小流域產(chǎn)流時空分布過程。
(1)平面混合產(chǎn)流模型。提出了山坡地貌水文響應(yīng)單元的劃分標準,將山洪地面響應(yīng)單元劃分為快速響應(yīng)單元(根據(jù)差異性又分為快、中、慢三種類型)、滯后響應(yīng)單元(根據(jù)差異性又分為快、中、慢三種類型),貢獻較小單元和其他單元等。研究不同山坡地貌響應(yīng)單元下墊面參數(shù)的差異,建立各類山坡地貌水文響應(yīng)單元與產(chǎn)流機制的對應(yīng)關(guān)系,實現(xiàn)超滲/蓄滿產(chǎn)流機制的平面上混合產(chǎn)流模型搭建。
(2)垂向混合產(chǎn)流模型結(jié)構(gòu)。時空變源混合產(chǎn)流模型的垂向結(jié)構(gòu)是綜合考慮了截留、填洼、超滲產(chǎn)流、蓄滿產(chǎn)流、優(yōu)先流、壤中流、基流和滲漏出流等過程。根據(jù)上述產(chǎn)流組分形成過程,基于非飽和下滲理論和水量平衡,采用不同的概念水庫來模擬土壤包氣帶和飽和帶,從上至下依次是表層土壤的毛細水庫、下層土壤的變動面積飽和產(chǎn)流水庫和地下水庫,各水庫產(chǎn)流的相關(guān)變量及參數(shù)見表1,垂向模型結(jié)構(gòu)如圖2。
(3)非飽和土壤下滲。非飽和土壤下滲采用Talbot和Ogden[20]提出的一種新的一維入滲和水分再分配的計算方法。將土壤含水量范圍離散成具有恒定含水量的區(qū)間,推導(dǎo)一維入滲方程以快速獲取濕潤鋒垂向位移,并考慮不同含水量單元之間的互吸力,實現(xiàn)再分配濕潤鋒橫向位移,實現(xiàn)高效穩(wěn)定的包氣帶土壤下滲計算方法,獲得包氣帶土壤的下滲過程;本文采用此方法進行小流域包氣帶土壤非線性下滲計算,該方法采用數(shù)值化方案,無需估算非線性梯度,使得最終的土壤下滲過程具有數(shù)值穩(wěn)定性。
表1 時空混合產(chǎn)流模型主要變量及參數(shù)
圖2 垂向混合產(chǎn)流模型結(jié)構(gòu)
2.2 FFMS水文模型軟件模塊化小流域暴雨洪水分析軟件FFMS由中國水利水電科學(xué)研究院劉昌軍博士負責(zé)研發(fā)。該軟件基于C++、FOTRAN和JAVA等語言開發(fā)完成,具有模塊化、參數(shù)化、智能化、可視化和自動化特點,軟件主界面見圖3。軟件集成了包括山洪災(zāi)害調(diào)查評價成果數(shù)據(jù)導(dǎo)入導(dǎo)出、小流域自動劃分及參數(shù)提取算法、模塊化建模環(huán)境、地貌水文響應(yīng)單元自動劃分、拖拽式手動和全自動建模、計算結(jié)果顯示等分布式水文模型的前后處理模塊,實現(xiàn)了自由組合建模和全自動模塊化建模功能;集成了多種國內(nèi)外分布式水文模型的氣象、蒸發(fā)、產(chǎn)匯流、坡面侵蝕、河道侵蝕和洪水演進等50余個算法模塊,具有參數(shù)自動率定和手動率定算法模塊等多個參數(shù)全局自動優(yōu)化算法。同時該軟件還集成了各種模型工具庫、機器學(xué)習(xí)參數(shù)庫和專家知識庫。詳細功能模塊見圖4。
圖3 FFMS軟件主界面
圖4 FFMS軟件架構(gòu)與功能模塊
3.1 控制方程本模型應(yīng)用動力波方法進行小流域坡面山洪過程的模擬計算。采用Godunov格式的有限體積法將計算區(qū)域進行空間離散,通過HLLC近似黎曼求解器計算單元界面上的質(zhì)量通量和動量通量[22]。本文應(yīng)用具有守恒格式的平面二維淺水方程見式(1)—式(2),鑒于模型和算法較為成熟,本文不再詳細介紹,模型詳細求解方法見文獻[22]。
式中:q為變量矢量,包括水深h,兩個方向的單寬流量qx和qy;g為重力加速度;u和v分別為x、y方向的流速;F和G分別為x、y方向的通量矢量;S為源項矢量;zb為河床底面高程;為謝才系數(shù),n為曼寧系數(shù)。
3.2 水動力學(xué)模型軟件簡介劉昌軍等采用OPENGL和C++語言開發(fā)了三維可視化水動力學(xué)模擬軟件FHMS,該軟件具有點云數(shù)據(jù)處理、DEM處理,網(wǎng)格剖分、參數(shù)設(shè)置、模型計算、參數(shù)率定和計算成果可視化展示等功能。本軟件可視化功能強,操作簡單,具有廣泛的推廣應(yīng)用前景,限于篇幅所限,不再詳細介紹,軟件主要功能界面見圖5。
4.1 流域基本情況以湖南寶蓋寺小流域為例,開展小流域暴雨山洪水文水動力學(xué)計算對比研究。小流域內(nèi)有寶蓋洞清水水文站,上游建有寒婆坳雨量站,小流域集水面積22.1 km2,河長9.1 km。小流域下墊面、土壤類型見圖6。選取2012年5月9日—5月10日的實測一場降雨徑流過程分別進行水文和水動力學(xué)模型對比計算分析。
圖5 三維水動學(xué)計算軟件主要界面
圖6 小流域土地利用和土壤質(zhì)地類型
4.2 水文建模與參數(shù)基于無人機獲取的1 m分辨率高精度DEM數(shù)據(jù),利用FFMS軟件全自動建立了寶蓋寺小流域分布式水文計算模型,共劃分15個模型計算單元,見圖7,模型參數(shù)采用2012年6、7月份的兩場實測暴雨洪水資料進行率定,得到此次模擬的模型參數(shù)取值,見表2。
表2 水文模型計算參數(shù)
圖7 小流域高精度DEM、計算單元劃分及坡面局部三維網(wǎng)格圖
4.3 水動力學(xué)建模與參數(shù)基于高分辨率DEM數(shù)據(jù),應(yīng)用FHMS軟件,建立了寶蓋寺小流域水動力學(xué)計算模型,計算網(wǎng)格模型采用非結(jié)構(gòu)化三角網(wǎng),共劃分網(wǎng)格單元1 764 215個,見圖7。該流域面積較小,地貌均質(zhì)性較好,以林地為主。糙率值是影響模擬結(jié)果的最重要參數(shù),以模型確定性系數(shù)為衡量標準,采用人工試錯法對糙率值進行率定,最終確定計算區(qū)域內(nèi)的坡面糙率值為0.09,溝道糙率值為0.03。
圖8 不同時間河道流量分布
圖9 不同時間的坡面淹沒水深
圖10 水文模型與水動力學(xué)方法計算值和水文站實測值
圖11 3個斷面水文模型與水動力學(xué)模擬洪峰流量過程
4.4 計算結(jié)果分析采用水文學(xué)和水動力學(xué)方法分別對該場次洪水進行了計算,水文學(xué)方法計算得到不同時間各河段洪峰流量平面分布見圖8,水動力學(xué)計算得到的河道淹沒水深見圖9。流域出口兩種方法的洪峰流量計算值和實測值見圖10。選取流域內(nèi)三個典型斷面對比分析兩種方法的洪峰流量見圖11,流域出口及各斷面計算的洪水峰現(xiàn)時間、洪峰流量結(jié)果見表3。
表3 洪水要素對比
(1)從圖8和圖9計算結(jié)果可以看出,水文模型和水動力學(xué)模型計算得到各子流域河道洪峰過程和洪峰流量基本一致。在降雨2 h后,9時左右在水文站下游河道開始形成徑流;10時小流域內(nèi)各支流河道均局部形成徑流,坡面局部地方也開始出現(xiàn)徑流;11時至12時河道內(nèi)洪水逐漸形成,并向下游演進。二維水動力學(xué)計算得到的坡面徑流產(chǎn)生、匯流、演進過程符合實際情況,與水文學(xué)方法計算結(jié)果一致。
(2)從圖10和表3可以看出,流域出口處兩種方法的模擬結(jié)果和水文站實測結(jié)果基本一致。此次洪水過程最大實測洪峰流量為27.7 m3/s,水文模型計算最大洪峰流量為31.6 m3/s,誤差約14%;水動力學(xué)模擬的最大洪峰流量為34.1 m3/s,誤差為23%。水文模型計算得到的峰現(xiàn)時間與實測值基本一致,水動力模型計算得到的峰現(xiàn)時間提前約1 h。
(3)從圖11中可以看出,三個典型斷面的水文學(xué)和水動力學(xué)計算結(jié)果基本一致,水動力學(xué)在3個斷面處均模擬出現(xiàn)了3次洪峰過程,而水文學(xué)模擬結(jié)果出現(xiàn)一次洪峰,過程曲線較為光滑,說明水動力學(xué)對洪峰過程的刻畫更加精細。
(4)本次采用的分布式水文模型在產(chǎn)流計算方面采用時空變源混合產(chǎn)流模型,此模型根據(jù)不同山坡地貌響應(yīng)單元的下墊面類型的差異,確定各類山坡地貌水文響應(yīng)單元所對應(yīng)的產(chǎn)流機制;建立并求解非飽和下滲曲線方程以快速獲取濕潤鋒垂向位移,考慮不同含水量單元之間的互吸力,實現(xiàn)濕潤鋒橫向位移再分配,得到包氣帶土壤水下滲過程;確定短歷時暴雨洪水的土壤最大含水量參數(shù),完成垂向和時段上超滲產(chǎn)流、淺層土壤蓄滿產(chǎn)流、深層土壤蓄滿產(chǎn)流機制時空轉(zhuǎn)換;有效解決了山丘區(qū)中小流域洪水預(yù)報面臨的水文數(shù)據(jù)精度低和產(chǎn)流機制時空變化問題,為無資料小流域的產(chǎn)流計算提供了新思路。
總的來看,水文學(xué)和水動力學(xué)方法對小流域暴雨山洪模擬結(jié)果都很好,說明自主研發(fā)的模塊化小流域洪水分析軟件和水動力計算方法均可用于暴雨山洪模擬。時空變源分布式水文模型模擬的結(jié)果與實際觀測結(jié)果更加接近;地表徑流二維水動力模型能夠處理實際流域的復(fù)雜地形,且可以處理具有間斷水流情況,具有較高的計算精度。
本文基于激光高精度激光點云數(shù)據(jù),應(yīng)用自主研發(fā)的時空變源混合產(chǎn)流模型、模塊化小流域分布式水文模型FFMS和二維水動力模型FHMS研究了寶蓋寺小流域坡面暴雨山洪過程。采用小流域出口水文站實測暴雨洪水資料驗證了水文學(xué)和水動力學(xué)方法計算結(jié)果,并對兩種方法計算暴雨山洪的效果、效率和精度進行對比分析,兩種方法均可用于小流域暴雨山洪計算分析,各有優(yōu)缺點,在山洪預(yù)報預(yù)警方面水文學(xué)方法更加簡單實用,在小流域暴雨洪水機理、演進過程和精細化分析方面水動力學(xué)更加精細和準確,本文提出的時空變源混合產(chǎn)流模型為無資料小流域暴雨洪水產(chǎn)流計算提供了新的研究思路,研發(fā)的水文學(xué)和水動力學(xué)計算軟件操作簡單、可視化功能強大、計算效率高,具有較大的推廣應(yīng)用價值。