MIKE21網格重編號的解決方法及應用

鄢忠清，熊海濱，胡小龍，楊緒海

（1.豐城市水利水電技術服務站，江西 豐城 331100；2.長江勘測規劃設計研究有限責任公司，湖北 武漢 430010；3.國家大壩安全工程技術研究中心，湖北 武漢 430010；4.流域水安全保障湖北省重點實驗室，湖北 武漢 430010；5.武漢大學水資源工程與調度全國重點實驗室，湖北 武漢 430072）

1 引言

MIKE21 軟件模型由丹麥水力學研究所（DHI）研制開發，模型整合了水動力模塊（HD）、對流擴散模塊（TR）、輸沙模塊（ST）、輸泥模塊（MT）、溢油模塊（OS）、粒子追蹤模塊（PT）等諸多模塊，廣泛應用于河流湖泊、河口海域的水流、波浪、泥沙及生態水質的模擬研究[1-6]。其中，水動力模塊作為核心模塊，是驅動其他模塊運行的基礎。前處理模塊中具有非常強大的自動網格生成器工具（Mesh generator），可對研究區域進行矩形網格、非結構網格劃分，并具有多種網格嵌套、局部加密等功能，對復雜邊界及水工結構物具有良好的適應性，各模塊的輸出面文件中每個網格單元都有特定的編號，且有對應的地理坐標。其中，在水動力模塊計算中，為提高計算速度，模型在默認條件下會對輸入地形文件中網格單元進行重新編號，導致計算輸出數據與前處理模塊輸出數據對應的網格編號不對應（相同編號網格對應不同坐標），且無明顯規律，從而使得計算模塊輸出結果無法直接與前處理模塊數據結合、交換。然而，基于MIKE21 研究河流工程領域眾多熱點問題需同時調用不同模塊的輸出數據，如水動力模擬中研究岸灘植被孳生導致河床阻力增加時需同時使用前處理模塊生成的地形高程數據和水動力模塊輸出的水深數據，又如河床沖淤分析需利用前處理模塊生成的地形高程數據及計算模塊輸出的網格面積數據[7]。因此，有必要對不同模塊間網格編號的差異性進行歸一化研究，使不同模塊數據能靈活交互，在充分利用MIKE 軟件各模塊強大功能的基礎上實現數據的有效整合。

2 方法介紹

基于地理坐標相同原則，提出了一種網格編號統一化的方法，克服了不同數據間無法直接有效結合使用的局限性，使不同網格編號數據能相互結合和交換，具體方法如下。

假定某套網格A 包含的網格單元數據集為Adata(xm,ym)，m∈[1,M]。其中，Adata(xm,ym)為第m個網格單元對應數據值；xm為第m個網格單元中心的橫坐標；ym為第m個網格單元中心的縱坐標；M為網格單元的數量。

另一套網格B 包含的網格單元數據集為Bdata(x′m,y′m)。其中，Bdata(x′m,y′m)為第m個網格單元對應數據值；x′m為第m個網格單元中心對應的橫坐標；y′m為第m個網格單元中心對應的縱坐標。

將2 套網格對應的網格編號及坐標歸一化：首先用升序提取網格A 數據集Adata(xm,ym)中的網格1～M編號的Adata(xm,ym)、xm、ym對應值，同理依次提取網格B 數據集Bdata(x′m,y′m)中網格1～M編號的Bdata(x′m,y′m)、x′m、y′m對應值。同時，將xm與x′m、ym與y′m從1～M依次進行匹配，共計2×M×M次，每當xm=x′r∩ym=y′r，r∈[1,M]時，使Adata(xm,ym)=Bdata(x′r,y′r)，可得與網格A 編號及坐標一致的新網格B 數據集Bdata(xm,ym)。

3 應用實例

植被孳生對河道水流條件影響是河流動力學領域的研究熱點之一。在進行水動力數值模擬時，糙率文件十分關鍵，模型率定、工況模擬均需要通過糙率的改變來實現。植被通常分布在地勢較高的岸灘上，且不同淹沒水深條件下植被對河道阻力的影響亦不相同。因此，在考慮植被阻力導致的河床糙率增加時，需同時結合前處理生成的高程數據及計算模塊輸出的水深數據，以長江中游河段某邊灘為例，基于本文提出的網格編號統一方法，結合地形及水文觀測資料，按照以下流程制作生成不同水深條件下的邊灘植被等效糙率文件，具體步驟如圖1所示。

圖1 糙率文件生成流程

（1）研究區域網格化和地形插值。包括研究區域邊界文件和地形文件的制作，采用前處理模塊的Mesh generator 工具首先導入邊界文件來對目標區域進行網格化，并將植被生長的岸灘范圍進行局部加密，設網格單元總數為N，然后導入特定年份的地形文件對網格單元進行地形插值，將插值后的地形分別導出成*.mesh 和*.dfsu 2 種文件格式，網格劃分如圖2所示。

圖2 研究區域網格劃分和地形插值

（2）岸灘范圍確定及特征值標記。首先復制1份步驟（1）所得*.dfsu文件，雙擊打開復制后的*.dfsu文件，利用工具欄下的多邊形選擇工具（Position of node selection polygon）選中步驟（1）中的加密區，在彈出的Edit Element Data 對話框中可以看到4 列數據，依次為網格單元編號、高程值、X坐標、Y坐標，數據行數為加密區的網格數量。然后復制高程1列數據到Excel軟件中，將其全部改為與整個河段高程范圍外任意一相同值，作為岸灘區標記，以便后續識別，最后將其復制到Edit Element Data 對話框中的高程1列替換原有高程值進行保存，如圖3所示。

圖3 岸灘加密區特征值標記

（3）網格單元的高程及坐標值提取。雙擊打開步驟（1）所得*.dfsu 文件，利用菜單欄Data 工具下的Select All 選中整個研究區域，在彈出的Edit Element Data 對話框中可以看到4 列數據，如圖4 所示，依次為網格單元編號、高程值、X坐標、Y坐標，其中單元編號為1～N，總計N行數據，將后3列數據復制到*.txt文件中備用；用同樣方法將步驟（2）中修改后的*.dfus文件的后3列數據復制到新的*.txt文件中。

圖4 網格地形高程數據（左、右圖分別截取了初始和末尾編號附近網格信息）

（4）特定流量下網格單元水深值的計算輸出。打開MIKE21 FM 模型的水動力模塊，在區域（Domain）導入步驟（1）的*.mesh 文件，設定模塊的進出口邊界，對模型參數進行率定和驗證，在輸出結果中選擇面文件（Area series）格式輸出，并在輸出項目（Output items）中勾選基礎變量（Basic variables）下的網格總水深（Total water depth），可得存儲網格水深的*.dfsu文件。

（5）網格單元水深及坐標值的提取。采用與步驟（3）相同的方式打開步驟（4）生成的*.dfsu 文件并選中整個研究區域，可以看到對話框中4 列含有網格編號、網格總水深、X坐標、Y坐標數據，如圖5 所示。對比圖4 發現，同樣的網格編號對應的坐標與步驟（3）中并不相同。將后3 列數據復制到h.txt 文件中備用。

圖5 網格水深數據（左、右圖分別截取了初始和末尾編號附近網格信息）

（6）網格高程數據及網格水深數據編號統一?；诒疚奶岢龅木W格編號統一方法，并通過Fortran編程技術實現，具體為將步驟（3）中生成的2個*.txt文件任意1個（2個文件網格編號統一）和步驟（5）生成的*.txt文件作為2 個輸入文件，建立實型數組Z（N）、X1（N）、Y1（N）和H（N）、X2（N）、Y2（N）分別讀取2個文件對應的3列數據。利用循環語句，首先將X1（1）與X2（1），X2（2），…，X2（N）比較，然后將X1（2）與X2（1），X2（2），…，X2（N）比較，依次類推，最后將X1（N）與X2（1），X2（2），…，X2（N）進行比較，共計N×N次比較。同理，對數組Y1（N）和Y2（N）進行比較，當X1（i）=X2（j）∩Y1（i）=Y2（j）時，其中i，j=1，2，…，N，將對應的H（j）對應的值賦給Z（i）。然后將運算后的結果數組Z（N）、X1（N）、Y1（N）輸出到h change.txt文件，即完成網格編號轉換，結果文件如圖6 所示，此時文件存儲的第一列數據為與步驟（3）中地形文件如圖3 所示的坐標對應的網格水深，且網格編號完全一致。

圖6 網格編號統一后的水深數據（左、右圖分別截取了初始和末尾編號附近網格信息）

（7）植被區等效河床糙率計算。根據植被分布區域和高程帶，再結合水深等植被生長特征參數確定河床等效糙率。包括將步驟（3）中生成的2個*.txt文件的高程列和步驟（6）輸出的*.txt 文件（h change.txt）的水深列復制到同一Execl 軟件工作表，使用其中的IF（IF（））函數進行雙重判定：首先基于步驟（2）的特征值判斷數據點是否在岸灘范圍內，其次根據高程值判斷該點的高程是否在植被生長帶，從而可以篩選出岸灘內植被生長帶的網格單元數據點。然后根據下式計算這些點的等效河床曼寧系數[8]。

式中：nv為植被區的等效曼寧阻力系數（s/m1/3）；k為二次流附加阻力系數；n為河床曼寧系數（模型率定值）（s/m1/3）；Cd為拖曳力系數；αv為形狀系數；hv為植被高度（m）；h為水深（m）；g為重力加速度（m/s2）；d為植株直徑（m）；cv為植被層植被體積與水體積之比[8]，計算公式為cv=πNd2/4；N為單位面積的植株數量（1/m2）；c為植被密度[8]。

將計算的等效糙率列數據復制替換到相同網格編號的步驟（1）或步驟（2）中的*.dfsu 文件中的高程值列即可實現結果可視化，文件可作為MIKE21 水動力模塊的糙率文件用于植被阻力數值模擬研究。不同流量下研究區域植被區等效糙率分布情況，如圖7 所示。由圖7 可以看出，與主河道相比，植被區所在河床糙率明顯增大；且不同位置水深的差異使得植被區內等效阻力空間上呈不均勻分布，最大糙率達0.205以上。由于在流量50000 m3/s時，植被已基本淹沒，隨著流量繼續增大，等效阻力反而減小，以上不同水深下植被區等效阻力的變化規律與類似研究一致，反映了糙率計算結果的合理性[9]。

圖7 不同流量下植被區等效阻力分布

4 結論

本文基于坐標相同原則，提出了解決MIKE21軟件不同模塊計算中網格編號存在差異的方法，克服了不同網格編號數據無法直接有效結合使用的局限性，并將其應用于植被等效糙率文件的生成，可有效利用MIKE21 軟件強大的網格處理功能，又能精確地反映不同水深下的植被阻力變化情況，從而為研究植被孳生對河道行洪、輸沙等功能影響研究提供基礎。