基于非均勻網(wǎng)格的高效高精度洪澇過程模擬

侯精明，王俊琿*，同 玉，張兆安，郭敏鵬，蘇 鋒，高徐軍

(1.西安理工大學(xué) 西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048；2.中國(guó)電建集團(tuán) 西北勘測(cè)設(shè)計(jì)研究院有限公司，陜西 西安 710065)

在全球氣候變暖背景下，極端暴雨發(fā)生的頻率和強(qiáng)度逐漸增加[1]，每年因極端暴雨洪澇災(zāi)害造成了巨大的人員傷亡和財(cái)產(chǎn)損失，嚴(yán)重制約經(jīng)濟(jì)社會(huì)的持續(xù)健康發(fā)展[2-3]。2016年汛期中國(guó)暴雨洪澇災(zāi)害頻發(fā)，長(zhǎng)江流域發(fā)生特大洪水，全國(guó)受災(zāi)人口達(dá)3 282萬，直接經(jīng)濟(jì)損失約1 470億元，與2000年以來同期均值相比，直接經(jīng)濟(jì)損失偏多51%[4]，水安全形式十分嚴(yán)峻。這就需要利用高精度數(shù)值模型研究洪水產(chǎn)生、演進(jìn)及消退過程，為政府決策部門提供行之有效的科學(xué)依據(jù)[5]。

數(shù)值模型在預(yù)測(cè)地表水流動(dòng)及相關(guān)過程中起著重要作用，其本質(zhì)是對(duì)非線性微分方程組進(jìn)行求解，其核心在于解的離散。解的離散程度會(huì)產(chǎn)生奇異性的解，奇異性使數(shù)值模擬受到阻礙。而網(wǎng)格的生成主導(dǎo)解的離散過程，因此一個(gè)高質(zhì)量的計(jì)算網(wǎng)格對(duì)模型的模擬精度和計(jì)算效率和有著很大的影響。目前，通常使用兩種網(wǎng)格對(duì)計(jì)算域進(jìn)行離散，即結(jié)構(gòu)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格。非結(jié)構(gòu)網(wǎng)格在處理復(fù)雜邊界問題時(shí)具有明顯的優(yōu)勢(shì)，但非結(jié)構(gòu)化網(wǎng)格關(guān)聯(lián)的數(shù)據(jù)結(jié)構(gòu)更復(fù)雜，控制方程的離散化更麻煩[6]。結(jié)構(gòu)化網(wǎng)格結(jié)構(gòu)簡(jiǎn)單，數(shù)值實(shí)現(xiàn)方便，且計(jì)算精度更高，已被眾多淺水水流模型采用[7-8]。

然而，在實(shí)際的洪水模擬過程中，計(jì)算域內(nèi)往往包含地形變化顯著的狹窄河道區(qū)域及面積廣闊的泛洪區(qū)[9]，且模擬通常涉及較大的研究區(qū)域。在均勻網(wǎng)格上，大范圍的高精度數(shù)值模擬不可避免會(huì)涉及大量計(jì)算節(jié)點(diǎn)，導(dǎo)致計(jì)算效率降低，而非均勻網(wǎng)格可以很好適用于具有復(fù)雜地形特征的區(qū)域。對(duì)地形變化明顯以及相對(duì)重要的區(qū)域提供精細(xì)網(wǎng)格，而其余區(qū)域則由粗網(wǎng)格離散，從而避免非必要的計(jì)算，減少計(jì)算成本。

近年來，隨著計(jì)算機(jī)技術(shù)的發(fā)展，網(wǎng)格劃分技術(shù)趨于成熟，周建中等[9]將自適應(yīng)網(wǎng)格技術(shù)應(yīng)用于潰壩洪水演進(jìn)模擬中，計(jì)算效率有一定提高，然而自適應(yīng)網(wǎng)格可能無法同時(shí)滿足質(zhì)量守恒和C-property[10]。質(zhì)量守恒影響水量，進(jìn)而影響水位。C-property通常反映模型保持所謂的良好平衡狀態(tài)的能力，即通量項(xiàng)和坡源項(xiàng)之間是平衡的，從而不會(huì)出現(xiàn)假動(dòng)量[6]；Liang[11]提出了非均勻網(wǎng)格方法以模擬洪水演進(jìn)過程，模擬精度高，但研究區(qū)域土地利用單一，且算例多為理想算例，不適應(yīng)于土地利用復(fù)雜（常為5種不同土地利用類型）、實(shí)際大規(guī)模范圍的城市發(fā)達(dá)地區(qū)洪澇模擬。基于此，本文提出一種新的網(wǎng)格劃分方法，分別對(duì)不同土地利用設(shè)置不同劃分準(zhǔn)則，以適應(yīng)實(shí)際大規(guī)模區(qū)域的洪澇模擬。此外，靜態(tài)非均勻網(wǎng)格可同時(shí)滿足質(zhì)量守恒與C-property[6]，對(duì)保證模型的數(shù)值穩(wěn)定性具有重要意義。

本文采用一種具有層次拓?fù)潢P(guān)系的、結(jié)構(gòu)化的靜態(tài)非均勻網(wǎng)格進(jìn)行淺層水流模擬的網(wǎng)格優(yōu)化，以高程梯度變化為判定變量調(diào)整計(jì)算域內(nèi)網(wǎng)格密度，在地勢(shì)平坦區(qū)域適當(dāng)加大網(wǎng)格尺寸，在水流模式受復(fù)雜地形特征影響的區(qū)域內(nèi)適當(dāng)加密網(wǎng)格并進(jìn)行高分辨率計(jì)算。與基于GPU加速技術(shù)的高分辨數(shù)值模型相結(jié)合，從而在保持解的精度基礎(chǔ)上，降低計(jì)算成本，提高模型的計(jì)算效率，以提升其在實(shí)際應(yīng)用中的性能。

1 非均勻網(wǎng)格劃分

本文采用遞歸層次的結(jié)構(gòu)網(wǎng)格（圖1），網(wǎng)格大小以網(wǎng)格劃分層級(jí)level反映，（記為lev），葉網(wǎng)格（如網(wǎng)格（i-1，j）和（i，j））劃分層級(jí)lev為0，其大小為Δx0和Δy0。對(duì)葉網(wǎng)格進(jìn)行細(xì)化可以得到子網(wǎng)格，子網(wǎng)格大小為Δx=Δx0/2lev，Δy=Δy0/2lev。為降低模型復(fù)雜性，網(wǎng)格與其鄰居網(wǎng)格大小滿足2∶1規(guī)則，不允許任何網(wǎng)格擁有大于或小于兩倍的鄰居，采用網(wǎng)格之間的相對(duì)關(guān)系來確定存儲(chǔ)鄰居網(wǎng)格信息[6]。

圖1 結(jié)構(gòu)化非均勻網(wǎng)格圖Fig.1 Structured non-uniform grid

生成精確反映區(qū)域地形特征的優(yōu)化的非均勻網(wǎng)格，一般遵循6個(gè)步驟：

1）定義一個(gè)初始均勻細(xì)網(wǎng)格來反映所考慮區(qū)域的地形數(shù)據(jù)的最高分辨率，網(wǎng)格大小為Δx，如圖2所示；在初始精細(xì)網(wǎng)格的基礎(chǔ)上設(shè)置虛擬背景粗網(wǎng)格，其大小為初始均勻細(xì)網(wǎng)格的整數(shù)倍，如圖3所示。一般粗化水平lev分為2級(jí)，最大粗化水平levm由式（1）計(jì)算：

圖2 根據(jù)地形生成的初始精細(xì)網(wǎng)格Fig.2 Initial fine grid generated from topographic data

圖3 定義levm=2的背景粗網(wǎng)格Fig.3 Defined background coarse grid with levm=2

式中，Δx、Δx0分別為初始網(wǎng)格和虛擬背景粗網(wǎng)格的單元大小。

2）對(duì)于任意背景粗網(wǎng)格，計(jì)算其內(nèi)部所有初始細(xì)網(wǎng)格的高程梯度（高程的1階導(dǎo)數(shù)）。對(duì)于初始細(xì)網(wǎng)格（i，j），其x和y方向的高程梯度與河床高程平均梯度分別為：

式中，zb為河床高程，Δx、Δy為初始細(xì)網(wǎng)格的大小，Grad_x、Grad_y分別為x和y方向的高程梯度，Grad為河床高程平均梯度。

3）對(duì)于任意背景粗網(wǎng)格，計(jì)算其內(nèi)部所有初始細(xì)網(wǎng)格的高程梯度變化（高程的2階導(dǎo)數(shù)），對(duì)于初始細(xì)網(wǎng)格（i,j），其x、y方向的高程梯度變化與河床高程平均梯度變化分別為：

式中，Gvar_x、Gvar_y分別為x和y方向的高程梯度變化，Gvar為河床高程平均梯度變化。

4）對(duì)于任意背景粗網(wǎng)格，取背景粗網(wǎng)格中所有初始精細(xì)網(wǎng)格的高程梯度變化最大值Gvar_max作為該背景粗網(wǎng)格的坡度變化值GvarΦ；取背景粗網(wǎng)格中所有初始精細(xì)網(wǎng)格對(duì)應(yīng)的土地利用類型數(shù)量最多者作為該背景粗網(wǎng)格的土地利用類型。

5）根據(jù)研究區(qū)域內(nèi)不同土地利用類型（即：建設(shè)用地、道路、裸地、林地和草地）重要程度設(shè)置不同閾值優(yōu)化網(wǎng)格（若僅有一種土地利用，則僅設(shè)一種閾值）。對(duì)于任一土地利用類型，設(shè)置兩層閾值grcr1和grcr2，即網(wǎng)格最高粗化水平為2級(jí)，當(dāng)背景粗網(wǎng)格的梯度變化值GvarΦ大于閾值grcr1，初始精細(xì)網(wǎng)格第一次粗化；當(dāng)背景粗網(wǎng)格的梯度變化值GvarΦ大于閾值grcr2，初始精細(xì)網(wǎng)格第二次粗化。否則，網(wǎng)格不需粗化。此外，可根據(jù)需求特定某區(qū)域，對(duì)該區(qū)域網(wǎng)格單獨(dú)細(xì)化。

6）引入2∶1規(guī)則保證網(wǎng)格過渡平穩(wěn)，即網(wǎng)格大小必須為鄰居網(wǎng)格的2倍、1倍或1/2倍。在此基礎(chǔ)上，定義網(wǎng)格順序編號(hào)、連接規(guī)則等[6]。最終，生成優(yōu)化后的非均勻笛卡爾網(wǎng)格，該網(wǎng)格鄰居信息完全由簡(jiǎn)單代數(shù)關(guān)系決定，避免使用顯式的數(shù)據(jù)結(jié)構(gòu)來存儲(chǔ)鄰居信息[12]，從而使存儲(chǔ)需求最小化。網(wǎng)格優(yōu)化情況如圖4所示。

圖4 優(yōu)化后的結(jié)構(gòu)化非均勻網(wǎng)格Fig.4 Optimized structured non-uniform grid

2 數(shù)值求解

2.1 控制方程

本文利用高精度數(shù)值模型模擬研究區(qū)域內(nèi)的洪澇災(zāi)害過程，模型所使用的控制方程為平面2維淺水方程，忽略運(yùn)動(dòng)黏性項(xiàng)、紊流黏性項(xiàng)、風(fēng)應(yīng)力和科氏力的2維非線性淺水方程矢量形式如下：

式中：x、y、t分別為笛卡爾空間坐標(biāo)與時(shí)間坐標(biāo)；h為水深；qx、qy為單寬流量；u、v表示x、y方向上流速；F、G表示x、y方向上通量矢量；S為源項(xiàng)矢量，包括降雨或下滲源項(xiàng)i、底坡源項(xiàng)及摩阻力源項(xiàng)；zb為河床底面高程；Cf為河床糙率系數(shù)，Cf=gn2/h1/3，n為曼寧系數(shù)。

2.2 數(shù)值求解

2.2.1 數(shù)值離散

采用中心格式的有限體積法，將變量存于網(wǎng)格中心，將式（4）在控制體 Ω上積分得：

式中，Г為控制體邊界，n為垂直于面元dГ的單位向量，F(xiàn)(q)n為相應(yīng)界面的通量。在笛卡爾坐標(biāo)系下，式（7）中的曲面積分項(xiàng)可以轉(zhuǎn)化為：

式中， Δx與 Δy分別為笛卡爾坐標(biāo)系下的網(wǎng)格單元（i，j）的左右和上下邊長(zhǎng)，F(xiàn)E、FW、FN和FS分別為網(wǎng)格單元（i，j）的東西南北4個(gè)界面的界面通量。

不同于均勻網(wǎng)格，非均勻網(wǎng)格的單元尺寸存在差異，如圖5所示。單元ic的東界面的通量由FE=FE1+FE2計(jì)算得出，其中FE1、FE2是通過兩個(gè)較小單元的西界面的中點(diǎn)的通量?？紤]以i為中心的較小單元，通量計(jì)算可能涉及位置w、i和e處的流量變量。由于守恒的流量變量h，qx和qy存儲(chǔ)在單元中心，因此，i處的流量變量可直接獲得，w和e處流量變量可通過自然鄰居內(nèi)插法獲得。

圖5 非均勻網(wǎng)格中通過控制體積邊界的流量Fig.5 Fluxes through the control volume boundary under the non-uniform grid

對(duì)界面通量進(jìn)行計(jì)算采用近似Godunov方法的HLLC近似Riemann求解[13]，該格式在處理干單元時(shí)的功能要優(yōu)于其他格式，能較好地處理復(fù)雜地形下干濕交替變化，具有較強(qiáng)的激波捕捉能力。其界面通量求解過程為：

式中：SL、SR、SM分別為計(jì)算單元左、右兩側(cè)及中間的波速，具體計(jì)算見式(10)。當(dāng)SL≥0和SR≤0時(shí)，計(jì)算單元界面的通量值分別由左、右兩側(cè)的水力要素確定；當(dāng)0≤SL≤SM和SM≤0≤SR時(shí)，計(jì)算單元界面的中間對(duì)流通量值由HLLC近似Riemann解給出，通量值詳細(xì)計(jì)算見參考文獻(xiàn)[13]。

式中，uL、uR分別為計(jì)算單元左、右兩側(cè)的流速，hL、hR分別為計(jì)算單元左、右兩側(cè)的水深，中間變量h*與u*由下式計(jì)算。

2.2.2 2階數(shù)值重構(gòu)

考慮到采用基于HLLC近似Riemann求解界面通量時(shí)在空間上僅具有1階精度，為使數(shù)值解的空間精度提高到2階，同時(shí)為了避免數(shù)值重構(gòu)后，在水面梯度大的地方出現(xiàn)數(shù)值振蕩的現(xiàn)象，模型采用TVDMUSCL方法進(jìn)行數(shù)值重構(gòu)[11]，表達(dá)式為：

其中：

式中，U為計(jì)算變量，φ為限制器函數(shù)，本模型采用Min mod限制器，其函數(shù)表達(dá)式為：

此外，為保證數(shù)值解整體上提高到2階精度，同時(shí)維持?jǐn)?shù)值解的穩(wěn)定性，對(duì)時(shí)間步長(zhǎng)采用兩階Runge-Kutta格式[14]來獲得2階精度。干濕動(dòng)邊界處理方面，水深容差為0.000 001 m以區(qū)別干濕網(wǎng)格單元[15]。底坡項(xiàng)使用模型作者提出的底坡通量法處理，摩阻項(xiàng)的計(jì)算使用分裂點(diǎn)隱式法以提高穩(wěn)定性[16]，計(jì)算過程中，Courant取0.5。

2.2.3 GPU加速方法

考慮到洪水過程模擬尺度較大，模型利用GPU并行計(jì)算技術(shù)，通過CUDA語言，將求解的淺水方程各項(xiàng)在GPU上運(yùn)行，在空間上實(shí)現(xiàn)大規(guī)模計(jì)算各網(wǎng)格水力要素信息；采用C++語言，將數(shù)據(jù)的讀寫、變量初始化在CPU上運(yùn)行。從而，在不降低精度條件下實(shí)現(xiàn)高速運(yùn)算[17]。圖6為GPU計(jì)算的流程圖。圖7為基于非均勻網(wǎng)格的數(shù)值模型進(jìn)行洪水模擬流程圖。

圖6 GPU計(jì)算流程圖Fig.6 Flowchart of the computation using GPU

圖7 基于非均勻網(wǎng)格的洪澇模擬流程圖Fig.7 Flowchart of flood simulation with non-uniform grids

3 非均勻網(wǎng)格在洪澇模擬中應(yīng)用

3.1 3個(gè)駝峰潰壩波模擬

本文選擇3個(gè)駝峰潰壩波模擬的理想算例來驗(yàn)證模型穩(wěn)定性。該潰壩問題被廣泛應(yīng)用于測(cè)試模型在2維復(fù)雜地形下處理非恒定流的能力[18]，具有一定的代表性。研究區(qū)域長(zhǎng)75 m，寬30 m，糙率n=0.018 s/m1/3，底部高程定義如下：

初始均勻網(wǎng)格分辨率為0.5 m，共計(jì)150×60個(gè)方形單元，以此為基礎(chǔ)進(jìn)行非均勻網(wǎng)格劃分，網(wǎng)格最高粗化水平為2級(jí)，設(shè)置兩層閾值：grcr1=0.001，grcr2=0.011。經(jīng)非均勻網(wǎng)格系統(tǒng)劃分，優(yōu)化后的非均勻網(wǎng)格數(shù)為3 216，研究區(qū)域地形及非均勻網(wǎng)格分布如圖8和9所示。

圖8 3個(gè)駝峰區(qū)域地形圖Fig.8 Topographic map of three humps

圖9 3個(gè)駝峰地形非均勻網(wǎng)格分布圖Fig.9 Non-uniform grid based on three hump terrain

算例共分兩部分：首先，驗(yàn)證和諧性，在整個(gè)計(jì)算域內(nèi)設(shè)初始水位為1.875 m，1 000 s后，計(jì)算域內(nèi)水位保持不變，表明模型具有良好的靜水和諧性。其次，修改初始條件，在x=16 m處建壩，大壩上游初始水位為1.875 m，下游干河床。t=0 s時(shí)，大壩瞬時(shí)全潰，模擬運(yùn)行300 s。此外，為驗(yàn)證水量平衡情況，地形四周設(shè)置為閉邊界，且無下滲情況，初始水量為1 125 m3，模擬300 s后，水量基本保持不變，因此，水量平衡狀況良好。圖10為P1、P2和P3這3個(gè)點(diǎn)位基于均勻網(wǎng)格與非均勻網(wǎng)格模擬的水位情況。

由圖10可看出，P1、P2和P3這3個(gè)測(cè)點(diǎn)基于非均勻網(wǎng)格模擬水位變化情況與基于均勻網(wǎng)格模擬水位情況吻合度較高，變化趨勢(shì)一致，均方根誤差RMSE分別為2.24%、2.23%、1.15%。因此，基于非均勻網(wǎng)格的數(shù)值模型具有較高的模擬精度。圖11為基于非均勻網(wǎng)格模擬的t=6 s和12 s的2維與3維水深圖。

圖10 P1、P2、P3測(cè)點(diǎn)的水位模擬情況Fig.10 Simulated water level at gauge P1, P2, P3

圖11 不同時(shí)刻的2維與3維水深圖Fig.11 3D and 2D water depth maps at different times

考慮到地形關(guān)于直線y=15 m對(duì)稱，理論上結(jié)果應(yīng)保持對(duì)稱性，圖11呈現(xiàn)的水深模擬結(jié)果具有較好的對(duì)稱性，符合地表水流運(yùn)動(dòng)規(guī)律。在潰壩水波演進(jìn)過程中，水波與主峰及小峰之間不斷接觸碰撞，t=12 s水流繞過主峰在下游形成干濕界面，最終在經(jīng)過多次碰撞及摩擦阻力的影響下，水流運(yùn)動(dòng)趨于平穩(wěn)。

本算例證明，模型有效捕獲由3個(gè)駝峰和區(qū)域邊界相互作用的激波引起的復(fù)雜流動(dòng)以及干濕界面及水位敏感區(qū)，且在干濕邊界附近未發(fā)生明顯的數(shù)值擴(kuò)散，具有良好的和諧性。此外，基于非均勻網(wǎng)格的模擬計(jì)算可以提高模型模擬效率，模擬采用PC機(jī)，搭載顯卡NVDIA RTX2070執(zhí)行，基于均勻網(wǎng)格的GPU加速的數(shù)值模型運(yùn)行26.95 s可完成300 s潰壩水流演進(jìn)過程，采用基于非均勻網(wǎng)格進(jìn)行計(jì)算，模型運(yùn)行時(shí)間為10.43 s，計(jì)算效率約為均勻網(wǎng)格的2.7倍。

3.2 Toce潰壩波模擬

本模型為Toce河物理模型，該區(qū)域潰壩波洪水演進(jìn)模擬[12]常用于驗(yàn)證水動(dòng)力學(xué)數(shù)值模型的模擬精度和數(shù)值穩(wěn)定性，具有一定典型性。該物理模型長(zhǎng)約50 m，寬約10 m，中部有一空水庫(kù)。初始均勻網(wǎng)格分辨率為0.05 m，共計(jì)約21萬個(gè)方形單元，以初始均勻網(wǎng)格為基礎(chǔ)進(jìn)行非均勻網(wǎng)格優(yōu)化，網(wǎng)格最高粗化水平為2級(jí)，設(shè)置兩層閾值grcr1 = 0.09，grcr2 = 0.12。經(jīng)非均勻網(wǎng)格系統(tǒng)劃分，優(yōu)化后的非均勻網(wǎng)格數(shù)為87 030，研究區(qū)域地形及非均勻網(wǎng)格分布如圖12和13所示。

圖12 Toce河區(qū)域地形圖Fig.12 Topographic map of Toce River

圖13 Toce地形非均勻網(wǎng)格分布圖Fig.13 Non-uniform grid based on Toce

利用數(shù)值模型分別基于非均勻網(wǎng)格及均勻網(wǎng)格對(duì)潰壩波洪水演進(jìn)過程進(jìn)行數(shù)值模擬，干濕水深的網(wǎng)格單元判定條件為0.000 001 m，Courant設(shè)置為0.5，輸入條件包括地形文件、入流文件等。左側(cè)設(shè)為入流邊界，右側(cè)為出流開邊界，其余設(shè)置為閉邊界[12]。圖14為該研究區(qū)域入流流量過程線。圖15為P4、P9和P21這3個(gè)點(diǎn)位基于數(shù)值模擬及實(shí)測(cè)水位情況。此外，為驗(yàn)證水量平衡情況，除左側(cè)入流邊界外，其余邊界設(shè)置為閉邊界，研究區(qū)域無下滲，入流過程總水量約為18 m3，模擬180 s后，研究區(qū)域總水量基本保持不變，因此，水量平衡狀況良好。

圖14 入流流量過程Fig.14 Inflow flow process data

圖15 P4、P9、P21測(cè)點(diǎn)基于模擬及實(shí)測(cè)水位情況Fig.15 Simulated and measured water level at gauge P4,P9, P21

由圖15可看出，P4、P9、P21這3個(gè)測(cè)點(diǎn)基于非均勻網(wǎng)格模擬水位變化與均勻網(wǎng)格模擬情況及實(shí)測(cè)水位變化吻合度較高。與實(shí)測(cè)值相比，非均勻網(wǎng)格模擬水位情況的納什效率系數(shù)分別為2.3%、1.0%和0.9%。因此，基于非均勻網(wǎng)格的數(shù)值模型可很好模擬潰壩波洪水演進(jìn)情況。此外，由于入流邊界在西側(cè)，潰壩水流的洪水演進(jìn)過程由實(shí)際地形的西側(cè)逐漸經(jīng)過P4、P9、P21測(cè)點(diǎn)，導(dǎo)致P4、P9、P21這3個(gè)測(cè)點(diǎn)實(shí)測(cè)水位在同一時(shí)刻數(shù)值不同，即峰值水位逐漸降低，峰現(xiàn)時(shí)間逐漸延后；而P4、P9這2測(cè)點(diǎn)地形較高且靠近入流邊界，P21測(cè)點(diǎn)地勢(shì)較平坦，因此，P4、P9測(cè)點(diǎn)水位與入流過程趨勢(shì)基本相同，呈現(xiàn)先增后減趨勢(shì)，而P21測(cè)點(diǎn)水位達(dá)到峰值后逐漸趨于平緩。圖16為不同時(shí)刻洪水演進(jìn)過程。

圖16 t=40 s和120 s時(shí)潰壩水流演進(jìn)過程分布Fig.16 Distribution of dam break flow evolution process at t=40 s and 120 s

此外，基于非均勻網(wǎng)格的模擬計(jì)算可以提高模型計(jì)算效率。模擬采用PC機(jī)，搭載顯卡NVDIA RTX2070執(zhí)行，基于均勻網(wǎng)格的GPU加速的數(shù)值模型需要耗時(shí)145 s才能完成洪水演進(jìn)過程，相同條件下，基于非均勻網(wǎng)格則僅需約為64 s，計(jì)算效率約為基于均勻網(wǎng)格模擬的2.3倍。

3.3 城市內(nèi)澇模擬

本文以陜西省西咸新區(qū)灃西新城為研究區(qū)域，降雨多集中在7至9月，且該區(qū)域土地利用類型多樣、建筑物林立、下墊面情況復(fù)雜[19]，選此研究區(qū)域模擬城市內(nèi)澇過程具有一定的代表性，圖17為研究區(qū)域地理位置概況圖。

圖17 研究區(qū)域地理位置概況圖Fig.17 Overview of the geographical location of the study area

利用機(jī)載激光雷達(dá)技術(shù)對(duì)研究區(qū)域地形進(jìn)行測(cè)算，得到分辨率為1 m共計(jì)312萬均勻網(wǎng)格單元的DEM數(shù)據(jù)，圖18為研究區(qū)域數(shù)字高程圖、正射影像圖。

圖18 研究區(qū)域數(shù)字高程及正射影像圖Fig.18 DEM date and orthophoto image map in the area

考慮到計(jì)算區(qū)域范圍較大，初始均勻網(wǎng)格較多，導(dǎo)致計(jì)算量大、耗時(shí)長(zhǎng)、成本較高。故以初始化均勻網(wǎng)格為基礎(chǔ)，采用優(yōu)化的非均勻網(wǎng)格系統(tǒng)對(duì)計(jì)算域內(nèi)不同的土地利用類型分別進(jìn)行劃分。網(wǎng)格最高粗化水平設(shè)置為2級(jí)，對(duì)于5種不同的土地利用類型（建設(shè)用地、道路、裸地、林地和草地）分別設(shè)置兩層不同的閾值grcr1與grcr2。

對(duì)于高程相對(duì)較低的草地、裸地、道路設(shè)置同一種閾值，本文首先將閾值設(shè)置為grcr1=grcr1草地=grcr1裸地=grcr1道路=0.50，grcr2=grcr2草地=grcr2裸地=grcr2道路=0.55，當(dāng)網(wǎng)格上的高程變化大于grcr1時(shí)，網(wǎng)格第1次粗化，當(dāng)網(wǎng)格上的高程變化大于grcr2時(shí)，網(wǎng)格第2次粗化，達(dá)到最高粗化水平（葉網(wǎng)格水平）；對(duì)于林地、建設(shè)用地設(shè)置另一種閾值，首先將閾值設(shè)置為grcr1=grcr1林地=grcr1建設(shè)用地=0.65，grcr2=grcr2林地=grcr2建設(shè)用地=0.70，當(dāng)網(wǎng)格上的高程變化大于grcr1時(shí)，網(wǎng)格粗化1次，當(dāng)網(wǎng)格上的高程變化大于grcr2時(shí)，網(wǎng)格第2次粗化，達(dá)到最高粗化水平。經(jīng)劃分，非均勻網(wǎng)格數(shù)約為130萬。

采用高分辨率數(shù)值模型分別基于初始均勻網(wǎng)格及該閾值設(shè)置條件下的非均勻計(jì)算網(wǎng)格對(duì)研究區(qū)域的內(nèi)澇過程進(jìn)行精細(xì)化模擬，輸入條件包括地形文件、降雨數(shù)據(jù)、土地利用文件等，四周為開邊界且無入流，計(jì)算過程中庫(kù)朗數(shù)取0.5。其中，地形文件分別為初始均勻方形網(wǎng)格的DEM數(shù)據(jù)及非均勻計(jì)算網(wǎng)格數(shù)據(jù)文件；降雨數(shù)據(jù)為2016年8月25日西咸新區(qū)云谷10號(hào)樓氣象站實(shí)測(cè)降雨，累計(jì)降雨量66 mm，降雨過程見圖19。

圖19 2016-08-25場(chǎng)次實(shí)測(cè)降雨Fig.19 2016-08-25 measured rainfall

根據(jù)研究區(qū)正射影像圖，采用最大似然分類法分別對(duì)構(gòu)建的初始均勻方形網(wǎng)格單元及非均勻網(wǎng)格單元?jiǎng)澐殖?種不同的土地利用（建設(shè)用地、道路、裸地、林地和草地），并制作成土地利用文件。不同土地利用的穩(wěn)定下滲率及曼寧系數(shù)見參考文獻(xiàn)[20]。

利用模型模擬計(jì)算域從開始降雨到7.5 h的積水過程。可根據(jù)模擬結(jié)果獲取到該區(qū)域不同時(shí)刻的內(nèi)澇點(diǎn)位置分布、積水面積、積水深度等內(nèi)澇情況。由模擬結(jié)果可知，t= 4 h時(shí)，該區(qū)域內(nèi)澇積水量為最大值，圖20為基于均勻網(wǎng)格與非均勻網(wǎng)格劃分的模擬結(jié)果積水量最大時(shí)分布情況。圖20中共標(biāo)記出內(nèi)澇積水量最大時(shí)的4處積水區(qū)（如圖20中橢圓圈中區(qū)域所示），其內(nèi)澇影響較為嚴(yán)重。

圖20 基于非均勻與均勻網(wǎng)格的積水點(diǎn)位模擬情況Fig.20 Simulated inundation based on the uniform grid and the non-uniform grid

模擬過程中基于非均勻網(wǎng)格的積水位置與均勻網(wǎng)格及實(shí)際情況的發(fā)生位置大致吻合（4處）。將基于計(jì)算域劃分的均勻網(wǎng)格的模擬積水程度與非均勻網(wǎng)格的模擬積水程度與實(shí)測(cè)數(shù)據(jù)[19]進(jìn)行對(duì)比，對(duì)比結(jié)果如表1所示。其中，實(shí)測(cè)數(shù)據(jù)由陜西省灃西新城管委會(huì)實(shí)地測(cè)量獲得。

由表1可知，基于結(jié)構(gòu)化非均勻網(wǎng)格的模擬與在分辨率為1 m的均勻網(wǎng)格上生成的模擬及實(shí)際水位情況相匹配。與均勻網(wǎng)格及實(shí)際情況相比，整個(gè)計(jì)算域內(nèi)4處積水位置的積水面積加權(quán)平均誤差分別為2.94%、8.95%，表明基于非均勻模擬精度與均勻網(wǎng)格模擬精度及實(shí)際情況相似。此外，本次模擬在模擬采用PC機(jī)，搭載專業(yè)計(jì)算顯卡NVDIA Tesla P100執(zhí)行，均勻計(jì)算網(wǎng)格模擬共耗時(shí)約13 300 s，非均勻網(wǎng)格模擬共耗時(shí)約5 500 s。速度約為原來的2.4倍?？梢?，基于結(jié)構(gòu)化非均勻網(wǎng)格的GPU加速的淺水流動(dòng)模型在模擬內(nèi)澇過程中可在保證準(zhǔn)確度的基礎(chǔ)上可有效提高模型的計(jì)算效率，縮短模擬時(shí)間。

表1 基于非均勻與均勻網(wǎng)格模擬水位及實(shí)際對(duì)比Tab. 1 Comparison of the field data and the simulated water level based on the non-uniform and uniform grid

4 結(jié) 論

計(jì)算域內(nèi)地形的概化精度對(duì)模擬結(jié)果有著很大影響，本文建立了基于非均勻網(wǎng)格的GPU加速的數(shù)值模型以求解淺水方程。并將新模型應(yīng)用到理想化地形的潰壩算例和實(shí)際大尺度地形的城市內(nèi)澇算例中，得到以下結(jié)論：

1）采用具有層次拓?fù)潢P(guān)系的非均勻網(wǎng)格劃分技術(shù)，以地形高程梯度變化為判定變量，根據(jù)研究區(qū)域不同土地利用類型制定不同劃分準(zhǔn)則以自動(dòng)調(diào)整網(wǎng)格密度，對(duì)地形變化明顯及水力要素敏感區(qū)域提供精細(xì)網(wǎng)格，而其余區(qū)域由粗網(wǎng)格離散，從而適應(yīng)實(shí)際洪水模擬中研究區(qū)域復(fù)雜地形特征的需求。

2）對(duì)高分辨率數(shù)值模型采用GPU并行計(jì)算技術(shù)，在不降低模擬精度的條件下，大幅提升模型計(jì)算效率。并與非均勻網(wǎng)格劃分模型相結(jié)合以模擬地表水流動(dòng)過程。由于相對(duì)重要的區(qū)域仍被較準(zhǔn)確描述，但網(wǎng)格數(shù)量大大減少，模擬速度在GPU加速的基礎(chǔ)上進(jìn)一步提高。因此，新模型可在保證模擬精度的同時(shí)，進(jìn)一步提高模型的計(jì)算效率。

3）將新模型應(yīng)用到實(shí)例中，在3個(gè)駝峰及Toce潰壩模擬中，模型可有效捕獲地形和邊界相互作用的激波引起的復(fù)雜流動(dòng)以及干濕界面及水位敏感區(qū)，驗(yàn)證了模型穩(wěn)定性，且模擬速度分別為基于均勻網(wǎng)格模擬的2.7與2.3倍；在灃西城市內(nèi)澇模擬中，基于非均勻網(wǎng)格模擬的積水位置與均勻網(wǎng)格模擬位置及實(shí)測(cè)情況大致吻合，積水面積加權(quán)平均誤差分別為2.94%、8.95%，速度約為均勻網(wǎng)格模擬的2.4倍，進(jìn)一步證明基于非均勻網(wǎng)格模型的高效性。

本研究提出的非均勻網(wǎng)格方法可描述區(qū)域的復(fù)雜地形特征，并與高分辨率數(shù)值模型結(jié)合以可靠地求解淺水方程，運(yùn)算速度在GPU加速基礎(chǔ)上進(jìn)一步提高。與基于均勻網(wǎng)格的模型相比，模型具有相似的求解精度，但計(jì)算成本較低，計(jì)算效率明顯提高。因此，在實(shí)際淺水流動(dòng)模擬中具有更好的潛力，可有效解決大規(guī)模洪水過程模擬遇到的精度低、效率慢等問題，適用于潰壩波洪水演進(jìn)過程及土地利用復(fù)雜多樣的發(fā)達(dá)地區(qū)的城市內(nèi)澇過程模擬。