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

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

1 引言

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

2 方法介紹

基于地理坐標(biāo)相同原則，提出了一種網(wǎng)格編號統(tǒng)一化的方法，克服了不同數(shù)據(jù)間無法直接有效結(jié)合使用的局限性，使不同網(wǎng)格編號數(shù)據(jù)能相互結(jié)合和交換，具體方法如下。

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

另一套網(wǎng)格B 包含的網(wǎng)格單元數(shù)據(jù)集為Bdata(x′m,y′m)。其中，Bdata(x′m,y′m)為第m個網(wǎng)格單元對應(yīng)數(shù)據(jù)值；x′m為第m個網(wǎng)格單元中心對應(yīng)的橫坐標(biāo)；y′m為第m個網(wǎng)格單元中心對應(yīng)的縱坐標(biāo)。

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

3 應(yīng)用實例

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

圖1 糙率文件生成流程

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

圖2 研究區(qū)域網(wǎng)格劃分和地形插值

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

圖3 岸灘加密區(qū)特征值標(biāo)記

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

圖4 網(wǎng)格地形高程數(shù)據(jù)（左、右圖分別截取了初始和末尾編號附近網(wǎng)格信息）

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

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

圖5 網(wǎng)格水深數(shù)據(jù)（左、右圖分別截取了初始和末尾編號附近網(wǎng)格信息）

（6）網(wǎng)格高程數(shù)據(jù)及網(wǎng)格水深數(shù)據(jù)編號統(tǒng)一。基于本文提出的網(wǎng)格編號統(tǒng)一方法，并通過Fortran編程技術(shù)實現(xiàn)，具體為將步驟（3）中生成的2個*.txt文件任意1個（2個文件網(wǎng)格編號統(tǒng)一）和步驟（5）生成的*.txt文件作為2 個輸入文件，建立實型數(shù)組Z（N）、X1（N）、Y1（N）和H（N）、X2（N）、Y2（N）分別讀取2個文件對應(yīng)的3列數(shù)據(jù)。利用循環(huán)語句，首先將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次比較。同理，對數(shù)組Y1（N）和Y2（N）進行比較，當(dāng)X1（i）=X2（j）∩Y1（i）=Y2（j）時，其中i，j=1，2，…，N，將對應(yīng)的H（j）對應(yīng)的值賦給Z（i）。然后將運算后的結(jié)果數(shù)組Z（N）、X1（N）、Y1（N）輸出到h change.txt文件，即完成網(wǎng)格編號轉(zhuǎn)換，結(jié)果文件如圖6 所示，此時文件存儲的第一列數(shù)據(jù)為與步驟（3）中地形文件如圖3 所示的坐標(biāo)對應(yīng)的網(wǎng)格水深，且網(wǎng)格編號完全一致。

圖6 網(wǎng)格編號統(tǒng)一后的水深數(shù)據(jù)（左、右圖分別截取了初始和末尾編號附近網(wǎng)格信息）

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

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

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

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

4 結(jié)論

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