基于水動力模型的泥石流風險分析

張弘強， 李文娟， 李大鳴， 唐智超， 湯稀博， 唐穎

(1.國網(wǎng)湖南省電力有限公司電力科學研究院， 長沙 410007； 2.高效清潔發(fā)電技術(shù)湖南省重點實驗室， 長沙 410007； 3.湖南省湘電試驗研究院有限公司， 長沙 410007； 4.天津大學水利工程仿真與安全國家重點實驗室， 天津 300072； 5.國網(wǎng)江西省電力有限公司萍鄉(xiāng)市湘東區(qū)供電分公司， 萍鄉(xiāng) 337016； 6.河海大學能源與電氣學院， 南京 211100； 7.長沙理工大學水利工程學院， 長沙 410114)

泥石流是一種破壞力強的常見災害，因此逐漸引起了中外學者的關(guān)注。對區(qū)域泥石流危險度進行評估，對泥石流風險進行分析，對于減輕泥石流災害的損失具有重要意義[1-2]。Mizuyama等[3]提出了用于計算由泥石流形成的溝床的比降變化的一維運動方程，以及用于計算泥石流影響范圍的二維運動方程。O’Brien等[4]將適用于黏性泥石流的Bing模型以及適用于稀性泥石流的Bagnold模型做了合并處理，提出了一種通用的泥石流運動模型，即稀釋塑性模型。De Joode等[5]構(gòu)建了一種增加了暴雨影響的泥石流災害模型，并將其與地理信息軟件相結(jié)合。Martinez等[6]提出了一種近似三維的數(shù)值模擬方法，用于模擬碎石較多的泥石流。Zou等[7]提出一種基于水文響應單元定量評價泥石流區(qū)域風險的方法。國內(nèi)學者從20世紀末開始對泥石流危險區(qū)域進行劃分，并得到了一系列相當具有價值的成果。譚炳炎[8]對泥石流溝進行了實地調(diào)查和統(tǒng)計，最先對泥石流溝的嚴重程度進行了數(shù)量化綜合評價。劉希林等[9]通過分析和調(diào)查，確定泥石流危險度的主要危險因子和次要危險因子，根據(jù)關(guān)聯(lián)次序確定它們的權(quán)重，以此提出了泥石流危險度的定量計算公式，極大地推動了泥石流量化評價的發(fā)展。李大鳴等[10]利用模糊賦值法計算各個泥石流危險因子的權(quán)重，建立了天津薊縣泥石流區(qū)域預報模型，確定了天津薊縣泥石流區(qū)域危險度。陳飛飛等[11]對泥石流危險評價問題進行探討，指出了目前評價方法中存在的問題，并探討了相應的解決方法。王高峰等[12]利用MATLAB多元非線性統(tǒng)計方法建立白龍江流域泥石流一次最大沖出量預測模型。李麗敏等[13]提出基于布谷鳥算法優(yōu)化后的輕量梯度提升機的泥石流預測。當前的泥石流評價主要集中在泥石流的靜態(tài)評價，在泥石流動態(tài)評價方面涉及較少[11]。

為構(gòu)建泥石流動態(tài)評價機制，現(xiàn)以水動力模型為研究依據(jù)，以降水、來流等自然條件為泥石流發(fā)生的推動因子，加之地形、地質(zhì)、植被、人口等影響因子，構(gòu)建泥石流風險分析模型。為準確判斷泥石流災害危險區(qū)域分布，為減輕泥石流損害提供指導性技術(shù)，現(xiàn)選取廣西5個市的行政區(qū)域作為研究區(qū)域，通過構(gòu)建水動力數(shù)學模型，以降水、地形、地質(zhì)、植被、人口等作為泥石流的影響因子，采用河網(wǎng)分級編碼方法概化河網(wǎng)，采用粗糙集理論計算權(quán)重，從而進行泥石流風險分析，并將計算模擬結(jié)果與廣西地質(zhì)災害分布統(tǒng)計和2011年全州縣泥石流災害事件進行對比，驗證模型的合理性和可靠性。

1 模型建立

1.1 控制方程

若假設(shè)在水深方向的動水壓強分布規(guī)律與靜水壓強相一致，則平面二維連續(xù)方程[14]為

(1)

式(1)中：t為時間；h為水深；u、v分別為x、y方向上的速度平均值；q為源匯項。

采用萊布尼茲公式當中的分解形式，并對其積分項作合并，從而得到二維淺水運動方程為

(2)

式(2)中：n為糙率；g為重力加速度。

1.2 控制方程離散求解

采用有限體積法對式(2)作離散，對于時間層面上的計算方式，采取水位與流量交錯的計算方式。對于控制方程的通用形式，可表示為

(3)

式(3)中：ρ為水的密度；φ為廣義變量；Γ為φ所對應的廣義擴散系數(shù)；S為φ所對應的廣義源項。

至于時間段Δt與控制體積ΔV，對控制方程的通用式(3)作積分可得

(4)

式中：div為求矢量的散度；grad為求矢量的梯度。

對于式(4)在空間與時間上的積分處理，集中表現(xiàn)在以下4個方面。

(1)瞬態(tài)項積分。

(5)

(2)源項積分。

(6)

式(6)中：SC對應節(jié)點C上廣義變量所對應的廣義源項的大??；SP對應節(jié)點P上廣義變量所對應的廣義源項的大?。沪誔為對應節(jié)點P上廣義變量的大小。

(3)對流項積分。

(7)

式(7)中：u為速度；(ρuφA)e、(ρuφA)w、(ρvφA)n和(ρvφA)s分別為ρuφA在所對應節(jié)點e、w、n和s的大??；(ρu)e、(ρu)w分別為ρu在所對應節(jié)點e、w的大小；(ρv)w、(ρv)s為ρv在所對應節(jié)點w、s的大小；Ae、Aw、An和As分別為A在所對應節(jié)點e、w、n和s的大??；φe、φw、φn和φs分別為φ在所對應節(jié)點e、w、n和s的大小。

(4)擴散項積分。

(8)

在對連續(xù)方程、動量方程等作離散求解時，只需要把相對應的物理量作為φ代入式(4)～式(8)，即可得到相對應的離散方程。對于連續(xù)微分方程，其對任何多邊形網(wǎng)格的顯式離散表達式為

(9)

式(9)中：Hi為i號網(wǎng)格單元的模化水深；T為所處時刻；Ai為i號網(wǎng)格單元所占面積大?。籐ik為i號網(wǎng)格、k號通道的長度；Qik為i號網(wǎng)格、k號通道的單寬流量；qi為i號網(wǎng)格的抽排水量。

1.3 邊界條件確定

對于模型邊界條件的確定，按照以下3個方面進行選取。

(1)研究區(qū)域邊界條件。以桂林、賀州、梧州、來賓、柳州5個市作為研究區(qū)域，以外輪廓線作為模型區(qū)域邊界條件。

(2)降水邊界條件。將開始時刻陸域網(wǎng)格內(nèi)的水深大小設(shè)為零，采用曲面差值將測量站所采集到的點雨量數(shù)據(jù)轉(zhuǎn)化為面雨量數(shù)據(jù)，以此構(gòu)成模擬計算泥石流的降水邊界條件。

(3)下墊面邊界條件。將廣西的地形、地質(zhì)、植被、人口分布等納入計算范圍，生成泥石流模型的下墊面邊界條件。

1.4 臨界條件確定

從泥石流的產(chǎn)生條件與影響因子出發(fā)，基于降水強度與入滲率的比較，分析泥石流的起動過程，即：①降水強度小于入滲率，雨水入滲，無積水；②淺層土體開始飽和，降水強度等于或略大于入滲率，積水導致淺層徑流，開始沖蝕松散土體；③徑流沖刷松散土體，地下滲流導致土體間的吸引力減小、抗剪強度減小，從而導致剪切破壞，加之水流的不斷作用，最終造成泥石流災害的發(fā)生。

根據(jù)團隊之前研究成果，給出泥石流發(fā)生的臨界條件[15-16]為

hc=[(ρsatgHcosθ+ρgH)(tanφ-tanθ)-ρsgHtanφ+τ1+τ2]/[ρmgcosθ(2tanθ-tanφ)+ρgtanφ]

(10)

式(10)中：hc為臨界徑流深度；H為破壞面入滲高度；ρsat為松散土體的飽和密度大小；ρm為松散固體物質(zhì)的密度大??；ρs為坡面松散固體物質(zhì)密度；ρ為水的密度；τ1為土體性質(zhì)的抗沖蝕能力的概化；τ2為植被覆蓋的抗沖蝕能力的概化；θ為斜坡坡度；φ為內(nèi)摩擦角。

1.5 泥石流數(shù)學模型構(gòu)建

圖1 廣西數(shù)字高程模型Fig.1 Digital elevation model of Guangxi

廣西境內(nèi)山地較多，地勢較為陡峭，高低落差較大，加之廣西地區(qū)雨量充沛，降水分布一般呈現(xiàn)出山區(qū)丘陵地帶多、平原河谷地帶少的規(guī)律。尤其在夏季，迎風坡面降水多、背風坡面降水少的特點，無疑為泥石流的發(fā)生創(chuàng)造了先決條件。從廣西數(shù)字高程模型圖(圖1)中提取整個研究區(qū)域的邊界，并將研究區(qū)域內(nèi)的35個市縣的地理信息數(shù)字化，得到研究區(qū)域的行政分區(qū)(圖2)。

圖2 模型行政區(qū)域劃分Fig.2 Administrative division of model

1.5.1 網(wǎng)格劃分

對于網(wǎng)格劃分，采用不規(guī)則、無結(jié)構(gòu)網(wǎng)格劃分技術(shù)，對研究區(qū)域進行自動剖分、逐步加密與分塊優(yōu)化。為了能夠確定研究區(qū)域占矩形范圍的大小，以邊長為0.01°的曲面正方形網(wǎng)格對矩形范圍內(nèi)的網(wǎng)格進行自動劃分，同時分別對單元、節(jié)點與通道進行編號。研究區(qū)域網(wǎng)格單元大小為2 km×2 km，模型共生成網(wǎng)格22 900個、節(jié)點23 768個、通道46 667個，局部示意如圖3所示。

1.5.2 地形處理

采用ArcGIS對廣西數(shù)字高程模型數(shù)據(jù)做預處理，從中提取地理高程數(shù)據(jù)，在單元節(jié)點位置與單元中心位置作曲面差值，得出研究區(qū)域范圍內(nèi)的地形數(shù)據(jù)。對數(shù)學模型作地形插值，利用了面積乘積公式：

(11)

圖3 局部單元-節(jié)點-通道關(guān)系Fig.3 Partial unit-node-channel relations

式(11)中：f為等待作插值的點；fi為已知的點；Si為以待插值點作頂點的三角形的面積大??；i∈Sj、i∈Sk分別表示Sj、Sk三角形包含頂點為fi的已知點。

1.5.3 河網(wǎng)處理

圖4 河網(wǎng)水系分布Fig.4 Distribution of river network

從基本數(shù)據(jù)中提取河網(wǎng)信息與水系分布，得到研究區(qū)域范圍內(nèi)的河網(wǎng)水系圖(圖4)。聯(lián)立控制方程以求解河道水動力信息的方法過于復雜，本文將河網(wǎng)進行數(shù)字化處理，基于二叉樹河網(wǎng)編碼方法進行改進，對河網(wǎng)進行編碼分級，為下一步計算做好先后順序上的準備。

其中，相關(guān)定義為：①河道是指沒有出流與匯流的河流；②岔口是指兩個及其以上河道的交匯點；③河段是指組成河道的基本單元；④過流斷面是指河段的銜接處。

編碼原則為：將河道岔口的入流端認定為上游，將河道岔口的出流端認定為下游。根據(jù)河道、岔口與河段的出現(xiàn)順序，分別進行編碼。若河道的上游岔口只有1條河道，則該河道稱為1級編碼河道，其上游無河流匯入，下游的河道匯入岔口。每一處岔口，其下游河道級別編碼需要等待該岔口上游河道的級別作完編碼之后，再繼續(xù)作編碼。按照這個規(guī)律往下游作編碼，到某個岔口若僅剩唯一河道，則視為該河網(wǎng)水系的出流邊界，編碼過程如圖5所示。根據(jù)該編碼原則，對研究區(qū)域內(nèi)的雙線河網(wǎng)作編碼，處理結(jié)果共有河道82115段、河道級別最大值已達41，模型局部河網(wǎng)編碼如圖6所示。

圖5 河網(wǎng)編碼過程Fig.5 Process of river network coding

2 模型計算

2.1 降水因子

降雨對于泥石流的產(chǎn)生，具有明顯的推動作用，一方面它是構(gòu)成泥石流的組分，另一方面它也是大量松散土體運動的載體，進而結(jié)合并產(chǎn)生泥石流。本文研究從國家氣象局共享數(shù)據(jù)庫中提取了2008年6月8—17日共10日的全國范圍內(nèi)的逐日網(wǎng)格降水量數(shù)據(jù)，通過自編程序，計算得出模型研究區(qū)域內(nèi)所需的115個節(jié)點的降水數(shù)據(jù)。之所以選擇這10日的降雨過程，是因為該降雨過程覆蓋了本文全部研究區(qū)域且雨量較大，可作為典型暴雨過程。

結(jié)合當?shù)氐湫捅┯攴峙淝€，將共計10日的日降水過程分配到240 h的時降水過程，以滿足數(shù)學模型的水動力計算要求。同時，采用泰森多邊形法，通過從115個雨量站所采集到的降雨量，推導求出計算范圍內(nèi)的面雨量，從而作進一步的產(chǎn)、匯流與輸水過程的數(shù)值計算，進而可以得出研究區(qū)域范圍內(nèi)網(wǎng)格表面與河網(wǎng)內(nèi)的積水深度最大值，如圖7所示。

2.2 地形因子

模型研究區(qū)域范圍內(nèi)的地形條件，毫無疑問是泥石流形成的主要誘發(fā)因子之一。本文將坡度作為評價泥石流危險度的指標，利用上文所述的模型地形處理方法，提取等高線、坡度等相關(guān)高程信息。之后，根據(jù)面積乘積公式，在網(wǎng)格單元節(jié)點處與中心處作高程曲面差值，從而得出研究范圍內(nèi)的地形，其三維地形如圖8所示，地形高程等值線如圖9所示。通過處理得到的地形高程等值線圖的走勢與廣西數(shù)字高程模型地形圖在大體上保持一致，所以對研究區(qū)域范圍內(nèi)地形的實際變化情況有較為真實的反映。

2.3 地質(zhì)因子

廣西當?shù)氐牡刭|(zhì)條件，能夠為泥石流的產(chǎn)生提供大量的松散土體。從數(shù)據(jù)資料中采集整個廣西地區(qū)的地質(zhì)分布情況，在數(shù)值計算中將研究區(qū)域內(nèi)的砂巖、礫巖、泥巖、頁巖等19種不同的地質(zhì)巖組進

行編碼，對這些巖組作數(shù)字化處理，從而得到研究范圍內(nèi)各個網(wǎng)格單元上的巖性屬性。局部巖組數(shù)字化結(jié)果如圖10所示，將其與來賓市地質(zhì)分布資料進行比較，說明該地質(zhì)編碼方法具有一定可行性。

圖6 河網(wǎng)編碼局部放大Fig.6 Partly enlarged and view of river network

圖7 模型計算區(qū)域網(wǎng)格表面與河網(wǎng)內(nèi)最大積水深度分布Fig.7 Maximum water depth distribution of mesh surface and river network in model calculation area

圖8 模型三維地形Fig.8 3D model terrain

圖9 模型地形高程等值線Fig.9 Elevation contours of model terrain

圖10 局部地質(zhì)數(shù)字化分布Fig.10 Partial distribution of geological digitization

2.4 植被因子

通過Google Earth軟件中的廣西衛(wèi)星地圖，獲取計算區(qū)域內(nèi)的下墊面植被條件，植被衛(wèi)星遙感圖如圖11所示。采用三原色基本原理并結(jié)合編程，將紅、藍、綠這三種色素的數(shù)值提出，并以綠色數(shù)值作為采集植被覆蓋率的依據(jù)。同時，將綠色數(shù)值劃分為不同的級別，分別代表林地、草地、田地等，經(jīng)過數(shù)值計算得出研究區(qū)域內(nèi)的植被覆蓋數(shù)字化信息，如圖12所示。比較圖11與圖12可知，研究區(qū)域范圍內(nèi)的植被覆蓋在經(jīng)過具體數(shù)據(jù)化之后，在大體上能夠與衛(wèi)星遙感數(shù)據(jù)中的植被覆蓋保持相同，表明所采集的植被覆蓋率情況是合理可信的。

圖11 植被衛(wèi)星遙感圖Fig.11 Remote sensing map of vegetation

圖12 植被數(shù)字化結(jié)果Fig.12 Digitized results of vegetation

2.5 人口情況

將研究區(qū)域范圍內(nèi)的人口分布作評價人類活動影響的依據(jù)。在數(shù)學模型的計算過程中，使用距離加權(quán)方法，并與控制面積加權(quán)方法作相結(jié)合進行使用，把所有行政區(qū)域范圍內(nèi)的總?cè)藬?shù)分撥至各個網(wǎng)格當中，人口分配公式為

(12)

2.6 影響因子權(quán)重計算

Pawlak[17]提出了一個關(guān)于粗糙集的理論，其特點在于結(jié)果完全取決于原始數(shù)據(jù)，不受主觀因素影響。根據(jù)粗糙集理論[18-19]，結(jié)合基于知識依賴性與信息熵兩種條件權(quán)重的優(yōu)點，結(jié)合多方面因素分析了條件屬性對研究區(qū)域中確定、不確定分類子集的影響，采用李大鳴等[20]提出的一種復合的綜合客觀權(quán)重計算方法。使用該方法對各因子權(quán)重進行計算，可以客觀定量地評定泥石流危險性等級，彌補了先前兩種方法的不足，確定各因子權(quán)重的方法可以表示為

wi=μwαi+(1-μ)wβi

(13)

式(13)中：wi為客觀權(quán)值；wαi為以代數(shù)觀為依據(jù)的權(quán)值；wβi為以信息觀為依據(jù)的權(quán)值；μ為各個因子。通過計算，得出各個影響因子的權(quán)重值，如表1所示。

在模型的數(shù)值計算當中，將水流流速、斜坡坡度、地質(zhì)巖性、人口分布情況與植被覆蓋情況分別當作評判泥石流危險度的降雨、地形、地質(zhì)、人類活動與植被覆蓋的指標。在對影響因子作數(shù)字化處理的基礎(chǔ)上，將各個影響因子的權(quán)值乘以其對應的數(shù)字化值，從而得到風險值[9]。通過對模型調(diào)節(jié)，得出一個較為合理的泥石流危險性分級，如表2所示。因此，典型暴雨過程下的泥石流災害分布，如圖13所示。

表1 影響因子權(quán)值Table 1 Impact factor parameter

表2 泥石流危險等級Table 2 Debris flow dangerous level

3 模型驗證

通過國家氣象局資源庫收集日降雨數(shù)據(jù)集，選擇2011年5月1—10日的降雨數(shù)據(jù)，計算得出積水深度以及各通道過流量并以其作為水動力推動因素，與斜坡坡度、地質(zhì)巖性、人口分布情況與植被覆蓋情況共同對泥石流危險度做出風險分析。通過泥石流數(shù)學模型計算得到該驗證時間段內(nèi)的泥石流災害危險區(qū)域分布，如圖14所示；通過長期統(tǒng)計結(jié)果所得到的廣西地質(zhì)災害分布，如圖15所示。

圖13 泥石流危險區(qū)分布計算結(jié)果Fig.13 Calculation results of debris flow hazard zone distribution

圖14 驗證時間段內(nèi)泥石流危險區(qū)分布計算結(jié)果Fig.14 Calculation results of debris flow hazard zone distribution in verification period

圖15 地質(zhì)災害分布統(tǒng)計結(jié)果Fig.15 Statistical results of geological hazard distribution

由圖13可得，典型暴雨過程下，廣西地區(qū)泥石流災害的高度危險區(qū)集中分布于資源縣、臨桂縣、陽朔縣、荔浦縣、柳江縣、柳州市等行政區(qū)域，泥石流災害的危險區(qū)集中分布于全州縣、龍勝縣、灌陽縣、融水縣、賀州市等行政區(qū)域。由圖15可得，實際統(tǒng)計結(jié)果中，廣西地區(qū)泥石流多發(fā)區(qū)集中分布于資源縣、臨桂縣、陽朔縣、柳江縣等行政區(qū)域。因此，通過對比圖13與圖15可知，實際已發(fā)生泥石流災害的地區(qū)，全部位于數(shù)學模型計算所得到的泥石流高度危險區(qū)范圍之內(nèi)，兩者泥石流危險區(qū)大體分布基本相同，說明計算結(jié)果基本可信。

另外，廣西全州縣咸水鄉(xiāng)洛江村廣坑糟屯于2011年5月9日4時20分左右發(fā)生了山體滑坡泥石流災害。由圖14可得，在5月1—10日這段驗證時間內(nèi)，整個研究區(qū)域只有全州縣計算得出了泥石流危險區(qū)。并且，受災地區(qū)恰好落在計算所得到的泥石流危險區(qū)與高度危險區(qū)最為密集的區(qū)域范圍內(nèi)(圖14右上角紅點、藍塊最密集處)，這與計算模擬結(jié)果吻合較好，這也進一步驗證了計算成果的可信度。

綜上所述，無論是長期地質(zhì)災害分布統(tǒng)計結(jié)果，還是特定時段泥石流災害事件，都能表明基于水動力學模型的泥石流災害風險分析具有一定的合理性及可靠性。

4 結(jié)論

以流速、水深作為水動力條件并采用了河網(wǎng)分級編碼方法，以坡度、巖性、植被覆蓋、人口分布分別作為地形、地質(zhì)、植被、人口等影響因子的評價指標，采用了一種以粗糙集理論為依據(jù)的、更加客觀的綜合權(quán)重計算法，從而構(gòu)建了泥石流風險分析數(shù)學模型。采用廣西地區(qū)地質(zhì)災害分布統(tǒng)計結(jié)果與2011年全州縣泥石流災害事件進行模型驗證，結(jié)果表明，該泥石流風險分析數(shù)學模型具有較好的可行性與可信度，在泥石流風險分析的應用當中是合理且可靠的。

針對泥石流破壞力強、規(guī)模大、活動頻繁等特點，通過水動力模型對泥石流風險進行分析，可以較準確地判斷泥石流災害危險區(qū)域分布,為減輕泥石流損害提供指導性技術(shù)。這對預防泥石流災害，減輕泥石流對人類的損害具有重要意義。