基于DEM的河道斷面構造改進方法及洪水演進精度評估

張文婷，李祉璇，張行南，劉永志

(1.河海大學水文水資源學院，南京 210098；2.水文水資源與水利工程科學國家重點實驗室，南京 210098；3.南京水利科學研究院水文水資源研究所，南京 210009)

洪水災害是自然災害之一，常導致嚴重的經濟損失和人員傷亡[1-2]。為了減輕洪水帶來的危害，需要不斷加深對洪水的研究。隨著計算機科學的快速發展，越來越多的水動力學數值模型用于模擬洪水過程，如MIKE 、HEC-RAS(hydrologic engineering center′s river analysis system)等[3-5]，為防洪減災措施的科學制定提供參考和依據。

水動力學數值模型可以得到與實際洪水演進過程較為吻合的模擬結果，需要河道斷面資料用于建模，但在無河道斷面地區和不方便實施河道斷面測量的地區其應用受到了限制[6-7]。河道斷面數據缺乏問題在許多地區都存在，盡管可以在建模前補充測量，但耗時長，花費昂貴[8]。針對這一問題，一些學者[9-11]開始探索從數字高程模型(digital elevation mode，DEM)中提取河道橫斷面數據應用到洪水模擬中，為缺少河道橫斷面數據的河段洪水模擬提供方法。在模擬中，研究者[12-16]發現從DEM中提取河道斷面數據存在對河道地形的描述不夠準確的問題，這將會影響洪水模擬準確性，因此,一些研究者嘗試通過使用河道斷面數據修正的方法對DEM數據做出調整，將河道實測斷面數據整合到DEM中，從而得出理想的模擬結果[17-18]。然而實際洪水模擬中許多河段的斷面數據不易獲得或測量成本過高，對此，一部分研究者[19]嘗試概化河道斷面，使概化的斷面與實測斷面保持一致的水位-流量關系；也有研究者探索通過結合水面以上地形數據推求水面以下的地形，例如俞茜等[20]采用三次樣條插值并合成完整斷面，得到的模型結果有較高可信度。但目前缺乏對于DEM提取河道斷面用于水動力學建模效果的系統研究和分析。

為此，本文在缺少河道斷面資料的地區以高精度DEM數據為支撐，通過構造斷面的方式快速獲取能夠應用于水動力學數值模型的河道橫斷面。分別采用由DEM直接提取河道橫斷面(以下簡稱DEM斷面)、實際勘測測量河道斷面(以下簡稱實測斷面)和改進的基于DEM構造的河道斷面(以下簡稱構造斷面)，進行河道水動力學數值模型建模，利用水動力學數值模型分別進行洪水模擬試驗，包括在不同頻率的洪水或不同模型結構(HEC-RAS 1D,HEC-RAS 1D-2D)中，分析不同河道斷面生成和獲取方法對模擬結果的影響，為缺少實測河道斷面資料河段的洪水模擬提供參考。

1 基于DEM的河道斷面提取方法和構造方法

1.1 DEM數據直接提取河道斷面

使用的地形數據由美國國家高程數據集(National elevation dataset，NED)以DEM的形式提供，其生成方式為通過航空激光掃描獲取激光雷達點云數據，經過數據處理得到本試驗采用的3 m×3 m分辨率的DEM柵格數據，是研究區范圍內可公開獲取的最高分辨率數據。DEM 被用作地形數據的輸入來源以創建河道橫斷面等，在RAS Mapper模塊中將橫斷面數字化，以轉換為 HEC-RAS 模型所需的河道幾何形狀，所有幾何數據都使用柵格表面的高程值進行編碼，作為DEM斷面的地形數據導入模型計算。

1.2 改進的基于DEM數據構造斷面

基于DEM的河道橫斷面構造方法[21]以Bhuyian提出的DEM校正方法為基礎，在河道橫斷面較為規則的河道中通過增大DEM斷面面積來增強河道的過流能力?；驹硎抢肈EM斷面高岸點指定的坡度確定各斷面左右岸的交點(圖1)，并假定其為各斷面最低點的近似位置。對于單渠道河道，邊坡的坡度從高岸點處的最大值逐漸減小至零。斷面最低點的近似位置確定之后，各橫斷面由通過斷面最低點的鉛直線分為左半部分和右半部分，每一部分面積相當于拋物線及其對稱軸圍成的面積，拋物線的垂直高度為均勻流深度y。此構造方法為斷面構造的Bhuyian方法，以下簡稱為Bhuyian法。

比對研究區實際勘測河道橫斷面情況，Bhuyian法構造的斷面形態較為接近實測斷面，但其面積存在偏小的情況，造成過流能力不足，需要在此基礎上增加構造斷面的過流面積，使得構造斷面與實測斷面的過流能力更為接近。本文提出一種改進的斷面構造方法，提出橫斷面最低點高程修正系數k，對橫斷面最低點高程進行修正，斷面形態仍由左半部分和右半部分的拋物線組成，最終構造橫斷面的方式見圖1。

圖1 橫斷面構造Fig.1 Schematic of section construction

構造的橫截面的幾何特性為

(1)

式中：B為河道斷面頂寬，m；B′為斷面寬度，m；y為均勻流水深，m；Zi為DEM表示的高岸坡點高程，m；Z′T為DEM最低點高程，m；k為橫斷面最低點高程修正系數。

斷面面積為

(2)

濕周為

P=∑Pi=Pl+Pr

(3)

水力半徑為

(4)

利用曼寧公式[式(5)]確定河道各斷面的均勻流水深，這是一個求得參考流量的均勻流水深精確解的迭代過程[21]。

曼寧公式為

(5)

式中：Q為參考流量,m3/s；n為糙率，其取值在《水工設計規范》等設計規范的基礎上結合實際河道類型及以往研究模擬中的相關取值經驗[22]確定；A為斷面過水面積，m2；R為水力半徑，m；S為河道坡度。

河道坡度由沿著河流中心線的DEM上可用的高程確定，確定糙率時要考慮到局部特征和地形，這一過程是為了實現估計的斷面最低點的縱向坡度和河道坡度值之間的最小差異。通過計算得出斷面最低點的高程為(Z′T-ky)。本文所提改進的斷面構造方法以Bhuyian法為基礎，利用橫斷面最低點高程修正系數k改進構造斷面的橫斷面面積，以下簡稱為改進法。

2 洪水演進模擬模型計算原理

用于河道水力計算的兩個模型為HEC-RAS 1D和HEC-RAS 1D-2D，均由美國陸軍工程兵團(USACE)的水文工程中心(Hydrologic Engineering Center)研發[23]。

HEC-RAS 1D中所用到的理論基礎是基于Saint-Venant方程組的水流連續方程和動量方程[24]

(6)

(7)

式中：AT為過流斷面面積，m2；t為過流時間，s；x為水平距離，m；Q為斷面流量，m3/s；ql為單位長度的旁側入流，m2/s；V為沿水流方向流速，m/s；g為重力加速度，m/s2；z為河道水位，m；Sf為摩阻比降，依據曼寧公式求得

(8)

式中：n為糙率；R為水力半徑，m。

在一維非恒定流模擬中，HEC-RAS模型對式(6)和(7)采用四點隱式差分格式離散求解方程。

HEC-RAS 1D-2D模型由進一步擴展一維模型得到。通過建立側向連接，將二維流動區域耦合到一維斷面。

二維連續方程為

(9)

式中：H為水面高程，m；h為水深，m；u和v為分別為x和y方向平均深度上的流速，m/s；q為源項，m/s，表示外部流入(如降水)。

二維動量方程中x方向為

(10)

y方向

(11)

式中：vt為渦流黏度系數；cf為摩擦因數；f為科里奧利參數。

側向結構上水流的流動由堰流方程或二維流動方程確定，計算側向結構上水流流動的標準堰流方程為

dQ=C(yws-yw)3/2dx

(12)

式中：yws為水面高程；yw為側向結構的高程；C為堰流系數。

3 實例應用

選取美國佛羅里達州邁阿密市內主河道邁阿密河作為研究案例，該河流包括一條干流和南北兩條支流。

3.1 研究區概況

邁阿密河發源于佛羅里達國家公園大沼澤地，流經邁阿密的市中心，長達8.9 km，河口位置為北緯25°46′14″，西經80°11′06″，研究區面積38.8 km2。

邁阿密河流經的邁阿密市位于佛羅里達大沼澤和比斯坎灣之間，屬于熱帶季風氣候，平均最高氣溫一般為29～38 ℃，年降水量1 710 mm，其中大部分降水發生在6月至10月中旬。此外，該區域常受颶風影響，如2004年弗朗西斯(Frances)颶風、2017年艾爾瑪(Irma)颶風。

3.2 數據

研究區域內流量和水位的觀測數據來自南佛羅里達州水管理區(SFWMD)環境監測數據庫，這些數據集可作為水動力學模擬的邊界條件輸入HEC-RAS模型。根據Anuar等[25]的研究，斷面布設間距為河道水面寬度的10～20倍時對水位過程模擬結果影響較小，考慮到邁阿密河大部分河道較為平直，研究區域內自上游向下游以500 m為間距布設河道橫斷面，布設情況、驗證站位置及邊界所處位置見圖2。

圖2 斷面布設情況、驗證站位置及邊界所處位置Fig.2 The distribution of cross-sections and the location of the gauging station and the boundaries

采用本文提出的改進的斷面構造法和DEM數據直接提取法生成研究區河道斷面。經過當地實際典型河道斷面與Bhuyian法概化的橫斷面對比分析，發現以Bhuyian法概化的橫斷面其過流面積雖然有所改善，但與實測斷面面積仍相差較大，通過改進法進一步擴大河道橫斷面面積，當最低點高程修正系數取值為1.1時，各典型斷面的橫斷面面積相對誤差均降至30%以下，縮小了約7%(表1)；若繼續增大最低點高程系數的取值，部分橫斷面的河底高程將低于實測斷面，綜合考慮橫斷面形態特征和過流能力，將最低點高程修正系數k取值為1.1。選取位于河道上游的斷面1、支流與干流交匯后的斷面2、位于河道下游的斷面3以及分別位于南北兩條支流上的斷面4和斷面5，比較實測斷面、Bhuyian法和改進法構造的斷面，如圖3所示，改進法構造的斷面在表示的河道測深細節方面優于Bhuyian法構造的，更為貼近實測斷面。

表1 橫斷面面積的相對誤差Tab.1 Relative error of cross-sectional area %

圖3 不同斷面形式對比Fig.3 Comparison of cross-sections

3.3 水動力學模擬參數設置

運行HEC-RAS 1D模型需要河道橫斷面的地形資料、邊界條件和糙率值。地形數據為研究區數字高程模型(DEM)，分辨率為3 m。模擬2004年颶風Frances期間(2004年9月1日00:00至2004年9月8日00:00)、2017年颶風Irma期間(2017年9月7日00:00至2017年9月14日00:00)和2014年6月17日至2014年6月24日的洪水演進過程，上邊界條件為干流和兩條支流的入流流量，下邊界條件為各颶風期間的潮位，初始條件為上邊界所在斷面流量初始值，模型計算步長為20 s，模擬結果輸出時間間隔為1 200 s。

在水力學模擬中，糙率是一個重要參數[26]，需要對河道和河漫灘的糙率分別賦值。根據《水工設計規范》等相關文獻資料，河道糙率為0.015～0.035，河漫灘糙率為0.075～0.150，結合研究區河道受到人為干預的實際情況，糙率參數的最終取值見表2。

表2 河道橫斷面糙率取值Tab.2 Parameter (Manning coefficient )

HEC-RAS 1D-2D模型由HEC-RAS 1D模型耦合水流的二維流動區域得到，輸入模型的數據在一維模型基礎上增加了洪泛區的地形數據和糙率值，水流的二維流動區域糙率取值依據美國地質調查局(USGS)的土地利用數據(NLCD)，對不同土地利用類型分區處理。研究區內不同類型下墊面的糙率取值見表3。

表3 不同下墊面類型糙率取值Tab.3 Roughness values of different underlying surface types

本研究主要目的是研究河道橫斷面數據對洪水演進模擬的影響，因此,同一種水動力模型結構中的模擬設置僅有河道橫斷面地形數據不同，其他參數均相同。

4 結果與分析

4.1 河道水位結果和流量結果對比分析

利用HEC-RAS 1D模型和HEC-RAS 1D-2D模型對Frances颶風期間(2004年9月1日至2004年9月8日)和Irma颶風期間(2017年9月7日至2017年9月14日)的洪水過程以及2014年6月17日至2014年6月24日的低水位過程進行水動力學模擬，根據驗證站處斷面實測水位過程，對實測斷面的HEC-RAS 1D和HEC-RAS 1D-2D模型進行模型驗證。驗證結果表明，HEC-RAS 1D模擬的3次水位過程均方根誤差分別為0.04、0.03和0.04 m，HEC-RAS 1D-2D模擬的3次水位過程均方根誤差分別為0.07、0.04和0.07 m，模擬結果的均方根誤差都不超過0.1 m，模型模擬結果的偏差較小，可以在該區域進行準確的洪水模擬。將構造斷面輸入HEC-RAS 1D模型和HEC-RAS 1D-2D模型，其模擬水位過程見圖4。

圖4 驗證站處采用構造斷面時模擬水位與觀測水位對比Fig.4 Simulated water level and observed water level with modified cross-section at the gauging station

圖4顯示構造斷面模擬出的水位過程與實測水位過程基本吻合，計算出模擬水位過程與實測水位過程的誤差(表4)，用以比較不同水動力模型結構和不同斷面形態模擬洪水水位過程的差別，計算結果表明，DEM斷面模擬出的水位數值和變化趨勢均與觀測值相差較大。在采用DEM斷面的4種模擬方案中，計算水位的平均誤差最高可達3.57 m，均方根誤差明顯高于采用構造斷面的模擬方案。

表4 模擬水位結果及分析Tab.4 Summary of simulations 單位：m

分析模擬水位誤差大小與實際觀測水位的關系(圖5)，發現在HEC-RAS 1D模型和HEC-RAS 1D-2D模型中，實測斷面和構造斷面的模擬誤差穩定在較低的數值，沒有隨水位而變化的趨勢，DEM斷面的模擬誤差隨著水位的升高呈現不斷減小的趨勢，當驗證站處水位達到1.61 m時，DEM斷面的計算誤差已經小于0.1 m，而在低水位條件下，DEM斷面的模型計算結果與觀測值相差較大。

圖5 水位模擬誤差與實測水位關系Fig.5 Relationship between simulation error and measured water elevation

研究區內的驗證站為水位站，缺少實測流量資料驗證構造斷面的準確性，但采用實測斷面的HEC-RAS 1D模型經過驗證可以模擬該研究區的洪水過程，故而將構造斷面產生的流量過程與實測斷面產生的流量過程進行比較，見圖6。圖6顯示，在不同的水位條件下，HEC-RAS 1D模型中構造斷面產生的流量過程與實測斷面產生的流量過程平均誤差為0.05 m3/s，HEC-RAS 1D-2D模型中兩者分別產生的流量過程變化趨勢及峰現時間比較接近，洪峰流量平均相對誤差約為10%?？紤]到偏大的模擬洪峰流量在實際應用中更為安全，偏差在可接受范圍。

圖6 驗證站處采用構造斷面與實測斷面時模擬流量對比Fig.6 Comparison of simulated flow between the modified cross-sections and measured cross-sections at the gauging station

4.2 淹沒范圍對比

洪水的最大淹沒面積由HEC-RAS 1D-2D模型中的RAS Mapper模塊計算得到，圖7為HEC-RAS 1D-2D 模型模擬的2017年Irma洪水過程中相應的淹沒面積，在3種斷面形式中，由構造斷面和實測斷面對應的HEC-RAS 1D-2D 模型模擬出的淹沒面積和變化趨勢在第72 h前比較接近，幾乎同時出現最大淹沒范圍，兩者的淹沒面積僅在退水過程中略有差異，相差最大的時刻出現在第96 h，淹沒面積相差約20%；DEM斷面由于缺少足夠的過水面積，且洪水期間上邊界的來水流量較大，致使河道中的水不斷地漫溢，洪泛區內洪水淹沒面積呈擴大趨勢，洪水的最大淹沒范圍產生于模擬計算結束的時刻。

分別使用實測斷面、構造斷面和DEM斷面構建的HEC-RAS 1D-2D模型計算出洪水淹沒范圍，圖8為這3種斷面數據模擬計算的最大淹沒面積對比。使用基于DEM斷面的模型計算出的最大淹沒范圍較基于實測斷面的模型產生的最大淹沒范圍偏大，與基于DEM斷面的模型相比，使用基于構造斷面的模型計算出的最大淹沒范圍與使用基于實測斷面的模型結果十分接近。

圖7 不同斷面形態下的洪水淹沒面積對比Fig.7 Comparison of flooded areas of different cross-sections

圖8 最大淹沒面積對比Fig.8 Comparison of maximum submerged area

圖9顯示了HEC-RAS 1D-2D 模擬的Irma洪水過程中實測斷面和構造斷面在不同淹沒深度時對應的淹沒面積對比，在不同的時刻，構造斷面模擬的不同淹沒深度對應的面積總體與實測斷面的模擬一致，約90%淹沒范圍水深小于1 m。

圖9 實測斷面和構造斷面在不同淹沒深度時對應的淹沒面積Fig.9 Area of the measured cross-sections and modified cross-sections at different submerged depths

在HEC-RAS 1D-2D模型的4個洪泛區總計4 729個二維流動網格中選取6個淹沒典型點，將構造斷面與實測斷面的模擬的結果進行比較。比較對象是Irma颶風期間6個典型點的淹沒過程，包括洪水到達時間、淹沒水深和洪水歷時，結果見圖10。

圖10 實測斷面和構造斷面的淹沒典型點水深Fig.10 Water depth of submerged typical points in measured cross-sections and modified cross-sections

圖10顯示，在二維洪泛區中構造斷面產生的淹沒過程和實測斷面接近，HEC-RAS 1D-2D采用構造斷面時能夠較好地反映最大淹沒水深，具有與采用實測斷面時幾乎相同的洪水到達時間，雖然構造斷面模擬的退水過程略快于實測斷面，但是在大多數淹沒典型點退水時長相差不超過4 h?？偟膩碚f，通過典型點的淹沒過程比較，在洪水到達時間、淹沒水深和洪水歷時方面構造斷面能夠達到與實測斷面接近的模擬效果，表明采用構造斷面能夠在該區域進行準確的洪水模擬。

5 結 論

以美國佛羅里達州的邁阿密河下游為例，分析河道斷面數據分別來源于DEM直接提取法、實地測量和本文改進的河道斷面構造法時，對模型HEC-RAS 1D，HEC-RAS 1D-2D模型模擬效果的影響，主要結論如下。

DEM數據在精確水力建模方面存在局限性，由于激光雷達掃描地形時忽略了水面以下的地形，DEM數據對河床的描述是不充分的，DEM斷面提供的過流能力不足，使用DEM斷面計算出的水位與實際觀測水位過程不相符合，且模擬出的淹沒深度、淹沒范圍與使用實測斷面的模擬結果存在明顯偏差。

采用本文提出的改進斷面構造的方式對DEM水面以下的數據進行構造矯正，在洪水演進模擬中能夠達到與使用實地測量河道橫斷面數據相近的效果，優于河道橫斷面僅從DEM數據直接提取法。通過對DEM資料中的河道橫斷面數據進行改進，可以節約實地勘測測量需要投入的時間和精力，從而增加DEM數據在廣泛實踐中的用處。

本研究中應用案例為平原河網地區，該地區受人類活動干預較大，橫斷面形狀相比于有明顯沖淤變化的河段更為規則，且采用的DEM數據分辨率較高。后續將對國內具備條件的流域進行擴展，并嘗試在山區型河道或由于沖淤導致橫斷面形狀不規則的河道區域，開展橫斷面獲取方法對水動力模型計算結果的準確性和適用性的進一步研究。