劉昌軍，文 磊，周 劍，趙懸濤，郭 良，魏永強(qiáng)

（1.中國(guó)水利水電科學(xué)研究院 減災(zāi)中心，北京 100038；2.河海大學(xué) 水文水資源學(xué)院，江蘇 南京 210098；3.中國(guó)科學(xué)院 西北生態(tài)環(huán)境資源研究院，甘肅 蘭州 730000；4.華北水利水電大學(xué) 水利學(xué)院，河南 鄭州 450045；5.湖南省水利水電科學(xué)研究院，湖南 長(zhǎng)沙 410007）

1 研究背景

我國(guó)山洪災(zāi)害點(diǎn)多面廣、突發(fā)性強(qiáng)、破壞力大，且多發(fā)生在偏遠(yuǎn)山區(qū)、交通不便、通訊不暢。山洪災(zāi)害是嚴(yán)重威脅群眾生命財(cái)產(chǎn)安全的主要災(zāi)種之一，日漸成為防洪減災(zāi)的短板，防御工作一直是我國(guó)防汛和防災(zāi)減災(zāi)工作中的難點(diǎn)和薄弱環(huán)節(jié)［1］。近年來(lái)無(wú)資料小流域暴雨山洪模擬研究得到快速發(fā)展［2］，目前常用的是采用水文學(xué)方法和水動(dòng)力學(xué)方法［3］。基于分布式水文模型的小流域暴雨洪水模擬得到快速發(fā)展和廣泛應(yīng)用［4-5］。黃燕等［6］開(kāi)展了水文模型在山洪模擬中的比較應(yīng)用研究，研究結(jié)果發(fā)現(xiàn)現(xiàn)有水文模型在模擬短歷時(shí)山洪效果一般。

隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，采用水動(dòng)力學(xué)方法進(jìn)行暴雨山洪過(guò)程的模擬是今后研究方向和熱點(diǎn)［7-8］。平面二維數(shù)值模型已成為洪水演進(jìn)模擬的主流方法，但間斷水流、干濕邊界處理和計(jì)算效率仍存在較多難點(diǎn)，是急需解決的關(guān)鍵問(wèn)題［9］。潘佳佳等［10］認(rèn)為基于運(yùn)動(dòng)波和擴(kuò)散波等簡(jiǎn)化水動(dòng)力學(xué)模型因忽略了流動(dòng)的力學(xué)因子，不能較好的模擬山洪的形成與演化機(jī)理。張培文等［11］采用有限元法研究了降雨條件下坡面徑流和入滲問(wèn)題的耦合模擬；王嬌進(jìn)行了坡面徑流—入滲耦合的解析解及無(wú)網(wǎng)格數(shù)值模擬研究［12］；基于Godunov格式的以求解黎曼近似解的淺水?dāng)?shù)值模擬技術(shù)已逐漸成熟，該方法能模擬光滑解，又能自動(dòng)適應(yīng)復(fù)雜流態(tài)變化，并得到廣泛應(yīng)用［13-14］。張大偉等［15］提出了新的基于Godunov格式的地表徑流二維水動(dòng)力模型，使干濕邊界問(wèn)題得到很好解決。張奕虹等［16］提出了新的水深重構(gòu)方法，解決了復(fù)雜地形上暴雨山洪數(shù)值模擬常面臨的計(jì)算水深小于零、非物理大流速等問(wèn)題。闞光遠(yuǎn)等［17-18］在GPU并行算法方面作了大量研究，發(fā)現(xiàn)使用GPU并行算法可大幅提高計(jì)算效率。GPU、CPU等并行計(jì)算技術(shù)已經(jīng)在水文和水動(dòng)力學(xué)技術(shù)中得到初步應(yīng)用，具有廣闊的應(yīng)用前景［19］。

針對(duì)暴雨山洪模擬計(jì)算問(wèn)題，以湖南瀏陽(yáng)市寶蓋寺小流域2012年降雨徑流過(guò)程為例，本文采用筆者提出的時(shí)空變?cè)椿旌袭a(chǎn)流模型和模塊化小流域暴雨山洪分布式水文模擬模型以及基于Godunov格式的二維坡面徑流數(shù)值模擬模型對(duì)小流域暴雨山洪進(jìn)行了對(duì)比研究，從建模過(guò)程、計(jì)算效率和應(yīng)用等方面，對(duì)比分析了兩種方法計(jì)算精度和計(jì)算效率及計(jì)算結(jié)果合理性，系統(tǒng)分析了坡面產(chǎn)匯流機(jī)制和微地形對(duì)降雨徑流過(guò)程影響，為無(wú)資料小流域暴雨洪水計(jì)算和預(yù)報(bào)預(yù)警提供技術(shù)借鑒。

2 水文模型及軟件簡(jiǎn)介

在山洪災(zāi)害調(diào)查評(píng)價(jià)成果和實(shí)時(shí)水文氣象大數(shù)據(jù)基礎(chǔ)上，筆者提出了以實(shí)時(shí)衛(wèi)星、雷達(dá)和臺(tái)站等多源降水融合數(shù)據(jù)和預(yù)報(bào)降雨數(shù)據(jù)為驅(qū)動(dòng)，以調(diào)查評(píng)價(jià)小流域拓?fù)潢P(guān)系及基礎(chǔ)屬性為基礎(chǔ)，以7類(lèi)水文單元（流域、河段、節(jié)點(diǎn)、分水、水源、洼地、水庫(kù)）水力關(guān)系為核心，以歷史暴雨洪水和率定參數(shù)為先驗(yàn)知識(shí)，以機(jī)器學(xué)習(xí)和并行計(jì)算為手段，以降雨蒸發(fā)、產(chǎn)匯流、河道演進(jìn)、水庫(kù)調(diào)蓄等水文過(guò)程為主線，基于人工智能和大數(shù)據(jù)技術(shù)，構(gòu)建以小流域?yàn)橛?jì)算單元，集成模型庫(kù)、人工智能算法庫(kù)、參數(shù)庫(kù)和先驗(yàn)知識(shí)庫(kù)為一體的模塊化分布式水文模型（FFMS）［20］，模型計(jì)算流程和結(jié)構(gòu)見(jiàn)圖1。

圖1 分布式水文模型流程和結(jié)構(gòu)

2.1 時(shí)空變?cè)椿旌袭a(chǎn)流模型圍繞山丘區(qū)中小流域地形地貌多樣、產(chǎn)流機(jī)制時(shí)空變化復(fù)雜等問(wèn)題，劃分了山丘區(qū)山坡地貌水文響應(yīng)單元，采用基于一維入滲理論的包氣帶土壤非線性下滲計(jì)算方法，計(jì)算濕潤(rùn)鋒下移過(guò)程及入滲量，提出了超滲、蓄滿時(shí)段轉(zhuǎn)換準(zhǔn)則，建立了在平面、垂向、時(shí)段上時(shí)空變?cè)椿旌袭a(chǎn)流模型［21］，以精確模擬山丘區(qū)中小流域產(chǎn)流時(shí)空分布過(guò)程。

（1）平面混合產(chǎn)流模型。提出了山坡地貌水文響應(yīng)單元的劃分標(biāo)準(zhǔn)，將山洪地面響應(yīng)單元?jiǎng)澐譃榭焖夙憫?yīng)單元（根據(jù)差異性又分為快、中、慢三種類(lèi)型）、滯后響應(yīng)單元（根據(jù)差異性又分為快、中、慢三種類(lèi)型），貢獻(xiàn)較小單元和其他單元等。研究不同山坡地貌響應(yīng)單元下墊面參數(shù)的差異，建立各類(lèi)山坡地貌水文響應(yīng)單元與產(chǎn)流機(jī)制的對(duì)應(yīng)關(guān)系，實(shí)現(xiàn)超滲/蓄滿產(chǎn)流機(jī)制的平面上混合產(chǎn)流模型搭建。

（2）垂向混合產(chǎn)流模型結(jié)構(gòu)。時(shí)空變?cè)椿旌袭a(chǎn)流模型的垂向結(jié)構(gòu)是綜合考慮了截留、填洼、超滲產(chǎn)流、蓄滿產(chǎn)流、優(yōu)先流、壤中流、基流和滲漏出流等過(guò)程。根據(jù)上述產(chǎn)流組分形成過(guò)程，基于非飽和下滲理論和水量平衡，采用不同的概念水庫(kù)來(lái)模擬土壤包氣帶和飽和帶，從上至下依次是表層土壤的毛細(xì)水庫(kù)、下層土壤的變動(dòng)面積飽和產(chǎn)流水庫(kù)和地下水庫(kù)，各水庫(kù)產(chǎn)流的相關(guān)變量及參數(shù)見(jiàn)表1，垂向模型結(jié)構(gòu)如圖2。

（3）非飽和土壤下滲。非飽和土壤下滲采用Talbot和Ogden［20］提出的一種新的一維入滲和水分再分配的計(jì)算方法。將土壤含水量范圍離散成具有恒定含水量的區(qū)間，推導(dǎo)一維入滲方程以快速獲取濕潤(rùn)鋒垂向位移，并考慮不同含水量單元之間的互吸力，實(shí)現(xiàn)再分配濕潤(rùn)鋒橫向位移，實(shí)現(xiàn)高效穩(wěn)定的包氣帶土壤下滲計(jì)算方法，獲得包氣帶土壤的下滲過(guò)程；本文采用此方法進(jìn)行小流域包氣帶土壤非線性下滲計(jì)算，該方法采用數(shù)值化方案，無(wú)需估算非線性梯度，使得最終的土壤下滲過(guò)程具有數(shù)值穩(wěn)定性。

表1 時(shí)空混合產(chǎn)流模型主要變量及參數(shù)

圖2 垂向混合產(chǎn)流模型結(jié)構(gòu)

2.2 FFMS水文模型軟件模塊化小流域暴雨洪水分析軟件FFMS由中國(guó)水利水電科學(xué)研究院劉昌軍博士負(fù)責(zé)研發(fā)。該軟件基于C++、FOTRAN和JAVA等語(yǔ)言開(kāi)發(fā)完成，具有模塊化、參數(shù)化、智能化、可視化和自動(dòng)化特點(diǎn)，軟件主界面見(jiàn)圖3。軟件集成了包括山洪災(zāi)害調(diào)查評(píng)價(jià)成果數(shù)據(jù)導(dǎo)入導(dǎo)出、小流域自動(dòng)劃分及參數(shù)提取算法、模塊化建模環(huán)境、地貌水文響應(yīng)單元自動(dòng)劃分、拖拽式手動(dòng)和全自動(dòng)建模、計(jì)算結(jié)果顯示等分布式水文模型的前后處理模塊，實(shí)現(xiàn)了自由組合建模和全自動(dòng)模塊化建模功能；集成了多種國(guó)內(nèi)外分布式水文模型的氣象、蒸發(fā)、產(chǎn)匯流、坡面侵蝕、河道侵蝕和洪水演進(jìn)等50余個(gè)算法模塊，具有參數(shù)自動(dòng)率定和手動(dòng)率定算法模塊等多個(gè)參數(shù)全局自動(dòng)優(yōu)化算法。同時(shí)該軟件還集成了各種模型工具庫(kù)、機(jī)器學(xué)習(xí)參數(shù)庫(kù)和專(zhuān)家知識(shí)庫(kù)。詳細(xì)功能模塊見(jiàn)圖4。

圖3 FFMS軟件主界面

圖4 FFMS軟件架構(gòu)與功能模塊

3 水動(dòng)力學(xué)模型簡(jiǎn)介

3.1 控制方程本模型應(yīng)用動(dòng)力波方法進(jìn)行小流域坡面山洪過(guò)程的模擬計(jì)算。采用Godunov格式的有限體積法將計(jì)算區(qū)域進(jìn)行空間離散，通過(guò)HLLC近似黎曼求解器計(jì)算單元界面上的質(zhì)量通量和動(dòng)量通量［22］。本文應(yīng)用具有守恒格式的平面二維淺水方程見(jiàn)式（1）—式（2），鑒于模型和算法較為成熟，本文不再詳細(xì)介紹，模型詳細(xì)求解方法見(jiàn)文獻(xiàn)［22］。

式中：q為變量矢量，包括水深h，兩個(gè)方向的單寬流量qx和qy；g為重力加速度；u和v分別為x、y方向的流速；F和G分別為x、y方向的通量矢量；S為源項(xiàng)矢量；zb為河床底面高程；為謝才系數(shù)，n為曼寧系數(shù)。

3.2 水動(dòng)力學(xué)模型軟件簡(jiǎn)介劉昌軍等采用OPENGL和C++語(yǔ)言開(kāi)發(fā)了三維可視化水動(dòng)力學(xué)模擬軟件FHMS，該軟件具有點(diǎn)云數(shù)據(jù)處理、DEM處理，網(wǎng)格剖分、參數(shù)設(shè)置、模型計(jì)算、參數(shù)率定和計(jì)算成果可視化展示等功能。本軟件可視化功能強(qiáng)，操作簡(jiǎn)單，具有廣泛的推廣應(yīng)用前景，限于篇幅所限，不再詳細(xì)介紹，軟件主要功能界面見(jiàn)圖5。

4 兩種方法對(duì)比分析

4.1 流域基本情況以湖南寶蓋寺小流域?yàn)槔?，開(kāi)展小流域暴雨山洪水文水動(dòng)力學(xué)計(jì)算對(duì)比研究。小流域內(nèi)有寶蓋洞清水水文站，上游建有寒婆坳雨量站，小流域集水面積22.1 km2，河長(zhǎng)9.1 km。小流域下墊面、土壤類(lèi)型見(jiàn)圖6。選取2012年5月9日—5月10日的實(shí)測(cè)一場(chǎng)降雨徑流過(guò)程分別進(jìn)行水文和水動(dòng)力學(xué)模型對(duì)比計(jì)算分析。

圖5 三維水動(dòng)學(xué)計(jì)算軟件主要界面

圖6 小流域土地利用和土壤質(zhì)地類(lèi)型

4.2 水文建模與參數(shù)基于無(wú)人機(jī)獲取的1 m分辨率高精度DEM數(shù)據(jù)，利用FFMS軟件全自動(dòng)建立了寶蓋寺小流域分布式水文計(jì)算模型，共劃分15個(gè)模型計(jì)算單元，見(jiàn)圖7，模型參數(shù)采用2012年6、7月份的兩場(chǎng)實(shí)測(cè)暴雨洪水資料進(jìn)行率定，得到此次模擬的模型參數(shù)取值，見(jiàn)表2。

表2 水文模型計(jì)算參數(shù)

圖7 小流域高精度DEM、計(jì)算單元?jiǎng)澐旨捌旅婢植咳S網(wǎng)格圖

4.3 水動(dòng)力學(xué)建模與參數(shù)基于高分辨率DEM數(shù)據(jù)，應(yīng)用FHMS軟件，建立了寶蓋寺小流域水動(dòng)力學(xué)計(jì)算模型，計(jì)算網(wǎng)格模型采用非結(jié)構(gòu)化三角網(wǎng)，共劃分網(wǎng)格單元1 764 215個(gè)，見(jiàn)圖7。該流域面積較小，地貌均質(zhì)性較好，以林地為主。糙率值是影響模擬結(jié)果的最重要參數(shù)，以模型確定性系數(shù)為衡量標(biāo)準(zhǔn)，采用人工試錯(cuò)法對(duì)糙率值進(jìn)行率定，最終確定計(jì)算區(qū)域內(nèi)的坡面糙率值為0.09，溝道糙率值為0.03。

圖8 不同時(shí)間河道流量分布

圖9 不同時(shí)間的坡面淹沒(méi)水深

圖10 水文模型與水動(dòng)力學(xué)方法計(jì)算值和水文站實(shí)測(cè)值

圖11 3個(gè)斷面水文模型與水動(dòng)力學(xué)模擬洪峰流量過(guò)程

4.4 計(jì)算結(jié)果分析采用水文學(xué)和水動(dòng)力學(xué)方法分別對(duì)該場(chǎng)次洪水進(jìn)行了計(jì)算，水文學(xué)方法計(jì)算得到不同時(shí)間各河段洪峰流量平面分布見(jiàn)圖8，水動(dòng)力學(xué)計(jì)算得到的河道淹沒(méi)水深見(jiàn)圖9。流域出口兩種方法的洪峰流量計(jì)算值和實(shí)測(cè)值見(jiàn)圖10。選取流域內(nèi)三個(gè)典型斷面對(duì)比分析兩種方法的洪峰流量見(jiàn)圖11，流域出口及各斷面計(jì)算的洪水峰現(xiàn)時(shí)間、洪峰流量結(jié)果見(jiàn)表3。

表3 洪水要素對(duì)比

（1）從圖8和圖9計(jì)算結(jié)果可以看出，水文模型和水動(dòng)力學(xué)模型計(jì)算得到各子流域河道洪峰過(guò)程和洪峰流量基本一致。在降雨2 h后，9時(shí)左右在水文站下游河道開(kāi)始形成徑流；10時(shí)小流域內(nèi)各支流河道均局部形成徑流，坡面局部地方也開(kāi)始出現(xiàn)徑流；11時(shí)至12時(shí)河道內(nèi)洪水逐漸形成，并向下游演進(jìn)。二維水動(dòng)力學(xué)計(jì)算得到的坡面徑流產(chǎn)生、匯流、演進(jìn)過(guò)程符合實(shí)際情況，與水文學(xué)方法計(jì)算結(jié)果一致。

（2）從圖10和表3可以看出，流域出口處兩種方法的模擬結(jié)果和水文站實(shí)測(cè)結(jié)果基本一致。此次洪水過(guò)程最大實(shí)測(cè)洪峰流量為27.7 m3/s，水文模型計(jì)算最大洪峰流量為31.6 m3/s，誤差約14%；水動(dòng)力學(xué)模擬的最大洪峰流量為34.1 m3/s，誤差為23%。水文模型計(jì)算得到的峰現(xiàn)時(shí)間與實(shí)測(cè)值基本一致，水動(dòng)力模型計(jì)算得到的峰現(xiàn)時(shí)間提前約1 h。

（3）從圖11中可以看出，三個(gè)典型斷面的水文學(xué)和水動(dòng)力學(xué)計(jì)算結(jié)果基本一致，水動(dòng)力學(xué)在3個(gè)斷面處均模擬出現(xiàn)了3次洪峰過(guò)程，而水文學(xué)模擬結(jié)果出現(xiàn)一次洪峰，過(guò)程曲線較為光滑，說(shuō)明水動(dòng)力學(xué)對(duì)洪峰過(guò)程的刻畫(huà)更加精細(xì)。

（4）本次采用的分布式水文模型在產(chǎn)流計(jì)算方面采用時(shí)空變?cè)椿旌袭a(chǎn)流模型，此模型根據(jù)不同山坡地貌響應(yīng)單元的下墊面類(lèi)型的差異，確定各類(lèi)山坡地貌水文響應(yīng)單元所對(duì)應(yīng)的產(chǎn)流機(jī)制；建立并求解非飽和下滲曲線方程以快速獲取濕潤(rùn)鋒垂向位移，考慮不同含水量單元之間的互吸力，實(shí)現(xiàn)濕潤(rùn)鋒橫向位移再分配，得到包氣帶土壤水下滲過(guò)程；確定短歷時(shí)暴雨洪水的土壤最大含水量參數(shù)，完成垂向和時(shí)段上超滲產(chǎn)流、淺層土壤蓄滿產(chǎn)流、深層土壤蓄滿產(chǎn)流機(jī)制時(shí)空轉(zhuǎn)換；有效解決了山丘區(qū)中小流域洪水預(yù)報(bào)面臨的水文數(shù)據(jù)精度低和產(chǎn)流機(jī)制時(shí)空變化問(wèn)題，為無(wú)資料小流域的產(chǎn)流計(jì)算提供了新思路。

總的來(lái)看，水文學(xué)和水動(dòng)力學(xué)方法對(duì)小流域暴雨山洪模擬結(jié)果都很好，說(shuō)明自主研發(fā)的模塊化小流域洪水分析軟件和水動(dòng)力計(jì)算方法均可用于暴雨山洪模擬。時(shí)空變?cè)捶植际剿哪Ｐ湍M的結(jié)果與實(shí)際觀測(cè)結(jié)果更加接近；地表徑流二維水動(dòng)力模型能夠處理實(shí)際流域的復(fù)雜地形，且可以處理具有間斷水流情況，具有較高的計(jì)算精度。

5 結(jié)論

本文基于激光高精度激光點(diǎn)云數(shù)據(jù)，應(yīng)用自主研發(fā)的時(shí)空變?cè)椿旌袭a(chǎn)流模型、模塊化小流域分布式水文模型FFMS和二維水動(dòng)力模型FHMS研究了寶蓋寺小流域坡面暴雨山洪過(guò)程。采用小流域出口水文站實(shí)測(cè)暴雨洪水資料驗(yàn)證了水文學(xué)和水動(dòng)力學(xué)方法計(jì)算結(jié)果，并對(duì)兩種方法計(jì)算暴雨山洪的效果、效率和精度進(jìn)行對(duì)比分析，兩種方法均可用于小流域暴雨山洪計(jì)算分析，各有優(yōu)缺點(diǎn)，在山洪預(yù)報(bào)預(yù)警方面水文學(xué)方法更加簡(jiǎn)單實(shí)用，在小流域暴雨洪水機(jī)理、演進(jìn)過(guò)程和精細(xì)化分析方面水動(dòng)力學(xué)更加精細(xì)和準(zhǔn)確，本文提出的時(shí)空變?cè)椿旌袭a(chǎn)流模型為無(wú)資料小流域暴雨洪水產(chǎn)流計(jì)算提供了新的研究思路，研發(fā)的水文學(xué)和水動(dòng)力學(xué)計(jì)算軟件操作簡(jiǎn)單、可視化功能強(qiáng)大、計(jì)算效率高，具有較大的推廣應(yīng)用價(jià)值。