基于RS和GIS的大南川水庫潰壩災害模擬方法研究

韓麗蓉，張 濤

(1.青海大學地質系，青海西寧810016;2.青藏高原北緣新生代資源環境重點實驗室，青海西寧810016;3.浙江省電力建設有限公司，浙江寧波315000)

一、引 言

我國已經建設的一部分水庫大壩存在著安全系數低的問題，大壩一旦潰決，勢必將對下游造成毀滅性的災難［1-2］。GIS技術運用于洪水災害的評估始于20世紀90年代，但是很多研究和技術方法都比較粗略，很少考慮利用其空間查詢分析功能提高計算分析的速度與準確性。劉仁義、劉南提出了有源淹沒和無源淹沒兩種模型;劉小生、黃玉生提出了利用“體積法”計算洪水淹沒面積［3］;周品、李勇等利用DEM數據結合“體積法”的概念公式，得到淹沒面積［4］;常燕卿、張福浩提出將洪水在特定河段內的DEM簡化為一斜平面，把淹沒范圍歸結為研究區域原始DEM被一個斜平面切割［5］。洪水演進模型主要是基于種子蔓延算法［1-7］，其不足之處在于必須編寫程序比較復雜。

二、洪水演進模型

本文探討了一種無需編程便可快速計算水庫潰壩后下游淹沒區的方法，適用于潰壩源頭與淹沒區落差較大地區。該方法結合本流域數字高程模型(DEM)，采用分段平面模擬方法。分段平面模擬法的思想是以本流域沿途鄉鎮作若干過水點特征橫斷面，將本流域劃分為若干個淹沒小段，把洪水在相鄰橫斷面內的DEM簡化為一個斜平面，并根據水動力演進模型模擬得到各斷面淹沒水深;然后，根據DEM解析出的山谷線高程和橫斷面的淹沒水深，計算得到河道兩端橫斷面的洪水高程;最后，通過將洪水淹沒面DEM與研究區域原始DEM疊加，并進行填挖方處理，得到洪水淹沒區情況。其中，填方區域表示洪水淹沒區，填方區與挖方區交線即為淹沒區邊界，將填挖方DEM轉換成面要素后刪除挖方部分，保留的填方部分為淹沒區范圍。

具體處理時對水庫及下游作如下簡化:設水流的滑動摩擦因數為μ;水庫壩高為x0，壩長為y0;水庫距下游目標研究淹沒區平距為D，高差為h。如圖1所示。

圖1 淹沒演進模型

假設水庫大壩由于地震等影響遭受劇烈破壞，造成1/2潰決，則可計算得出潰口面積大小

由于壩址處被已潰壩體阻擋，水流不會瞬間落至1/2壩體處，假設水流瞬間從壩頂作自由落體至3/4壩高，水庫內水的原始速度為0，根據勢能轉換關系得到3/4壩高處的水流速度

由以上兩式可得水庫潰壩后庫容放空時間

式中，K為水庫的庫容量;L為洪水潰口處流量;M0為洪水潰口面積。

水流出庫后將沿著地表順流而下，地表近似為一傾斜面，則水流在地表流動的加速度為

式中，θ為斜面坡度;μ為地面對水流的滑動摩擦因數。

設水流到達淹沒區的時間為t，則根據運動方程有

經過計算得到水流到達淹沒區特征過水點橫斷面的時刻t與流速vt

考慮水庫潰壩后洪峰經過時沿途淹沒高度，設在t時刻洪水淹沒區的橫截面面積為Mt=xtyt;洪水流量恒定不變，為Lt=Mtvt;潰口處流量為L0=Lt;即有xtytvt=M0v0，從而得到t時刻洪峰到達的淹沒區的淹沒高度

三、大南川水庫潰壩災害分析

1.建立研究區域DEM

大南川水庫建于1974年，位于青海省湟中縣，集雨區面積407 km2，水庫庫容1300萬m3，壩頂高程 2 640.5 m，壩長460 m，最大壩高 46.5 m，屬于湟水支流南川河水系。利用研究區域的地形圖，在ArcMap中線性內插出高程數據，建立三角網TIN;然后把TIN格式轉換成GRID格式DEM。

2.潰壩模擬

(1)提取研究區域山谷線

從DEM提取山谷線，山谷線與斷面交點為該斷面最低點，進而根據山谷線，并以洪水高度為依據確定淹沒水位。首先利用ArcMap的柵格計算器得到與淹沒前DEM地形相反的反地形;再按山脊線的提取方法得到山谷線［8］。

(2)提取特征斷面

本文選擇大南川水庫流域的8個人口密集區為代表，即大南川河水庫上游1個及下游7個橫斷面，繪制橫斷面圖。斷面方向與南川河法線(山谷線的法線)基本重合，考慮到大南川水庫的有限庫容量及洪水流經的極限流量，經反復試驗，提取的斷面最高點與最低點的高差基本保持在45 m，保證了在洪水深度的計算中能獲得較高圖解精度的斷面長度。

本文選擇繪制7個分段縱斷面圖。(3)計算洪水到達時間與洪水水深由大壩數據及洪水演進模型分析可知，潰口大小為

此時水流速度為

水流流量大小為

水庫放空時間為

假設洪水的過水橫斷面近似為矩形，地面對水流的滑動摩擦因數μ=0.015，重力加速度g=10 m/s2，利用各鄉鎮橫斷面與縱斷面，并結合淹沒模型，可得到各橫斷面長度yt、各縱斷面坡度及坡角、洪水到達各橫斷面的時間、洪水到達時的流速及洪水深度等。

(4)計算淹沒后洪水DEM

根據以上計算得出的洪水深度，在ArcMap中模擬繪制出洪水淹沒面DEM。由于兩橫斷面間距較大，因此采用插值方法加繪洪水淹沒面等高線，等高線的走向與山谷線的法線基本重合;然后利用ArcMap的3D分析工具獲取等高線與山谷線交點處的高程屬性，并將交點處的山谷線高程值加上洪水深度，即得到該等高線所在位置的洪水面高程;最后把得到的洪水面高程賦予繪制的等高線，生成TIN格式的DEM，并將其轉換為分辨率為5 m的GRID格式的DEM。因為研究區域落差較大，最高點與最低點之差達400 m，因此如果以某特定高程面為洪水淹沒面會帶來很大的誤差，甚至不能得出結果，導致模擬失敗。

(5)確定淹沒區范圍

將上面得到的洪水DEM與研究區域原始DEM疊加，并進行填挖方處理，得到填挖方范圍。其中填方區域表示洪水淹沒區，填方區與挖方區交線即為淹沒區邊界。挖方區在本研究中無意義，應將其刪除。但由于在ArcGIS中不易對DEM進行裁剪處理，因此需要先將填挖方DEM轉換成面要素，然后刪除挖方部分，保留填方部分，即得到洪水淹沒區，如圖2所示。

圖2 洪水淹沒區范圍

(6)淹沒區域三維可視化

在ArcScene中以前面討論得到的研究區域DEM的值為高程值，在其表面疊加顯示遙感影像和淹沒區范圍［9］。由于ArcScene中無法對行政區地名進行標注，因此利用ArcScene的3D編輯器制作三維地名及洪水到達時間的三維符號，得到淹沒區的三維景觀，如圖3所示。

圖3 淹沒區域三維場景

3.淹沒區損失分析

(1)淹沒區遙感影像處理

首先，利用ArcGIS的數據格式轉換功能將面要素格式的淹沒區范圍轉換成柵格格式(像素大小為5 m)，再轉換為Erdas能用的有地理參考的IMG格式(其存在背景噪聲);其次，在ERDAS視窗中對淹沒范圍進行邊緣檢測，去除背景噪聲，不關閉視窗，在Vector菜單中選擇相應功能將淹沒范圍轉換成矢量格式，再轉換為柵格格式;再次，在ERDAS中以淹沒區域為裁剪區域，對研究區域遙感影像進行掩模裁剪，得到淹沒區范圍的遙感影像;最后，對淹沒區遙感影像進行非監督分類［1］，將淹沒區地物分成5類:房屋、耕地、林地、水面、空閑地，進而得到各類的面積。

(2)損失情況分析統計

淹沒區域主要覆蓋西寧市城南新區和城中區，面積合計151 000 000 m2，人口合計29.7萬。本次洪水淹沒面積為30 726 900 m2，估計受災人口為30 726 900÷151 000 000×29.7=6.0萬人。

四、總結與展望

本文綜合運用ERDSA和ArcGIS等手段，以西寧市湟中縣的大南川水庫為例，討論了洪水演進模型，提出了一種快速簡便的適用于落差較大地區洪水淹沒分析的洪水演進模型。通過該模型可以計算出下游鄉鎮遭遇洪水的時間、洪水淹沒的深度、洪水的流速等重要洪水參數，進而能確定洪水淹沒范圍。模擬演示了流域的三維場景，并對淹沒范圍的遙感影像進行非監督分類，估計得出潰壩洪水給下游造成的損失，效果比較理想。

本研究中還存在不少需要改進之處:①所建立的洪水演進模型只考慮了洪峰演進情況，沒有考慮庫容大小、上游來水和集雨區降水量，會給研究結果造成一定的誤差;②所使用的圖件成圖時間較早，分辨率偏低，影響研究結果的準確性。洪水災害的預測是一個非常復雜的課題，其受到很多因素的影響，還有待進一步研究。

［1］王韶玉.基于DEM的壩堤潰決洪水淹沒評價模型與方法研究［D］.武漢:華中科技大學，2010:56-58.

［2］周遠方.大南川水庫潰壩的數值模擬研究［D］.長沙:長沙理工大學，2010:63-64.

［3］劉小生，黃玉生.基于Arc/Info的洪水淹沒面積的計算方法［J］.測繪通報，2003(6):46-48.

［4］劉仁義，劉南.基于GIS的復雜地形洪水淹沒區計算方法［J］.地理學報，2001，56(1):1-6.

［5］葛小平，許有鵬，張琪，等.GIS支持下的洪水淹沒范圍模擬［J］.水科學進展，2002，13(4):456-460.

［6］何宗宜，黃春雄，韓用順.水庫下游洪水淹沒災害評估系統的研究［J］.測繪通報，2003(2):24-27.

［7］馮平，紀恩福，盧永蘭.水庫下游洪災淹沒損失的估算［J］.災害學，1995，10(1):5-9.

［8］湯國安，楊昕.地理信息系統空間分析實驗教程［M］.北京:科學出版社，2006:85-86.

［9］張成才，劉丹丹，于欣，等.基于ArcGIS的東平湖洪水淹沒場景三維可視化［J］.鄭州大學學報:工學版，2008，29(1):88-90.

［10］黨安榮，賈海峰，陳曉峰，等.遙感圖像處理教程［M］.北京:清華大學出版社，2010:187-192.