毛思媛，賈艷紅，假冬冬，鄭祥民，范寶山，張佳欣

(1.華東師范大學地理科學學院地理信息科學教育部重點實驗室，上海 200241；2.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210029； 3.中水東北勘測設計研究有限責任公司，吉林 長春 130061；4.水利部寒區(qū)工程技術研究中心，吉林 長春 130061)

隨著全球極端氣候事件的增加，洪水頻率也相應增加[1]，這不僅對生命財產(chǎn)安全構成了重大威脅，也對河岸生態(tài)系統(tǒng)產(chǎn)生了重大影響。Berthelot等[2]研究表明，洪水頻率的增加將導致不同洪水重現(xiàn)期內(nèi)物種的組成、多樣性及群落結構出現(xiàn)顯著差異，特別是在洪水頻發(fā)區(qū)，植被的密度、更新率與物種多樣性會顯著減少。為有效減輕洪澇災害造成的巨大損失，促進河流生態(tài)建設，模擬流域洪水演進并分析洪水淹沒對河岸生物群落的影響對恢復河流基本功能、改善河流水生態(tài)環(huán)境及促進河流生態(tài)治理具有重要意義。

近年來，隨著洪水管理及生態(tài)河道治理政策的推行，如長江流域“共抓大保護，不搞大開發(fā)，生態(tài)優(yōu)先，綠色發(fā)展”戰(zhàn)略、黃河流域“積極實施生態(tài)保護和高質(zhì)量發(fā)展”戰(zhàn)略，松花江流域的洪水管理與生態(tài)環(huán)境保護也日漸受到人們的重視。目前已在嫩江洪水管理與洪水演進方面展開了一些研究，如孫亞楠等[3]利用MIKE11水動力學方法模擬了不同洪水組合下嫩江受第二松花江洪水頂托的影響,楊浩淼[4]基于HydroInfo二維水動力模型模擬了嫩江灣濕地擬建道路對河道防洪的影響,李洋[5]采用HEC-RAS一二維水動力耦合模型對嫩江左岸齊齊哈爾以上防洪保護區(qū)進行了洪水風險圖模型的構建與驗證。在洪水對河岸生態(tài)系統(tǒng)的影響方面國內(nèi)也有部分研究成果，如饒品增等[6]利用Landsat 7/8數(shù)據(jù)和MNDWI方法對2013年黑龍江省洪水對植被的影響進行了評估；余曉等[7]通過MIKE11一維非恒定流水動力學模型模擬了額爾古納河洪水淹沒范圍并分析了濕地植被變化情況。目前關于不同頻率洪水運動特征及其對生態(tài)系統(tǒng)影響的研究相對較少。

不同重現(xiàn)期的洪水淹沒程度顯著影響著河岸植被的群落分布。本文在利用水動力學方法建立嫩江流域河道洪水演進模型的基礎上，進一步模擬了典型河段4種不同重現(xiàn)期的洪水運動特征，并選取MODIS衛(wèi)星的MOD13Q1植被指數(shù)數(shù)據(jù)說明不同重現(xiàn)期洪水運動特征及淹沒范圍具有的生態(tài)學意義，旨在一定程度上彌補當前河道管理與生態(tài)河道建設缺乏規(guī)范與定量性研究的不足，同時也為河道實現(xiàn)近自然治理的目標提供基礎數(shù)據(jù)。

1 研究區(qū)概況

嫩江是松花江的最大支流，位于黑龍江省西北部，地處溫帶大陸性季風氣候區(qū)，夏季高溫多雨，降雨集中，上游河網(wǎng)密度大，支流多，河道全長1 370 km，流域面積29.7萬 km2，是我國水旱災害頻發(fā)的地區(qū)[8]。嫩江江橋—大賚河段位于嫩江下游三岔河口以上平原地帶(圖1)，河長約225 km[9]，主要流經(jīng)黑龍江省齊齊哈爾市、大慶市與吉林省白城市，主河道布設有江橋、大賚兩個水文測站，以及位于兩水文測站間的白沙灘水位站(123°53′E, 46°18′N)。該河段為典型的平原性河道，地形復雜，主流寬度小，河灘寬度大，平均坡降約0.0045‰～0.005‰[9]。河道水流平緩，河槽蜿蜒曲折，河道內(nèi)沙洲、汊河密布，多呈網(wǎng)狀，且沿江兩岸灘地上廣泛分布著大片濕地、泡沼群及牛軛湖。大洪水期間，江道與泡沼相連，兩岸灘地洪水漫溢，最寬處可達20.5 km，最大水深達7.4 m，河槽調(diào)蓄量大，水位漲落緩慢，洪水歷時長，使得汛期兩岸堤防的防洪任務較重[10]，也顯著影響著沿岸地區(qū)的植被生長。

圖1 研究區(qū)概況

2 一、二維水動力耦合模型構建

MIKE FLOOD一、二維耦合模型是把一維模型和二維模型連接在一起進行動態(tài)耦合的模型系統(tǒng)。耦合模型既利用了一維模型和二維模型的優(yōu)點，又避免了采用單一模型遇到的網(wǎng)格精度和準確性問題[11]?；诖?，本文采用丹麥水資源及水環(huán)境研究所(DHI)研發(fā)的MIKE FLOOD進行河道洪水演進模擬，模型以非恒定流運動的圣維南方程組(Saint-Venant)作為河道洪水模擬計算的基本方程[12]。

2.1 MIKE11一維水動力模型的構建

MIKE11一維水動力模型[13]的建立主要由河網(wǎng)數(shù)據(jù)、斷面數(shù)據(jù)、邊界條件及水動力學參數(shù)等組成，具體參數(shù)設置如下：

a.河網(wǎng)數(shù)據(jù)：根據(jù)2003年繪制的嫩江流域1∶5 000地形圖，本文對江橋—大賚河段河流走向線進行描繪，并將其概化生成河網(wǎng)文件。

b.斷面數(shù)據(jù)：利用ArcGIS對高精度的DEM數(shù)據(jù)進行斷面提取，由于嫩江江橋—大賚河段河長超200 km，考慮到模型運行時的穩(wěn)定性及準確性，本文共選取65處典型斷面，各斷面間距約1 000～4 000 m，并對各斷面最低點與左右岸進行標注。

c.邊界條件：選取江橋站1998年、2013年與不同設計頻率(1%、2%、5%、20%)的洪水流量過程作為模型的入流邊界條件，選取大賚站同期的水位過程及水位-流量關系作為模型的出流邊界條件。由于研究河段長、洪水傳播時間長，本文以各年份汛期(6—10月)作為計算時段，歷時約123 d，各測站均采用日水位-流量資料。

d.水動力學參數(shù)：根據(jù)《松花江流域防洪規(guī)劃》和實地查勘情況，本文將初始水位設為130.65 m，曼寧系數(shù)初始設置為0.035，模擬時間步長為5 s，輸出結果保存為.res11文件。

2.2 MIKE21 二維水動力模型的構建

MIKE21 FM 水動力模塊的建立主要包含網(wǎng)格剖分、地形構建、糙率分區(qū)及水動力參數(shù)(干濕邊界、渦黏系數(shù))的設定等。

a.網(wǎng)格剖分：本文采用非結構三角網(wǎng)格對建模區(qū)域進行劃分，并根據(jù)研究區(qū)域地形起伏靈活確定網(wǎng)格大小，如對河堤處進行網(wǎng)格加密。本次模型共劃分網(wǎng)格96 886個，節(jié)點52 573個。

b.地形構建：將高程散點數(shù)據(jù)導入網(wǎng)格中進行插值，生成地形文件。

c.糙率分區(qū)：嫩江流域共有耕地、林地、草地、濕地、水體、人工表面及裸地7類地物，地物類別數(shù)據(jù)主要來源于全國地理信息資源目錄服務系統(tǒng)GlobeLand30數(shù)據(jù)集 (http://www.webmap.cn/commres.do?method=globeIndex)。糙率的選取主要參考《松花江流域防洪規(guī)劃》與陳曉更等[14]的研究，并結合現(xiàn)場實地查勘與工程經(jīng)驗，在模型率定過程中進行相應的調(diào)整。

d.水動力參數(shù): 由于計算區(qū)域中存在隨水位漲落而變化的動邊界，為保證模型計算的連續(xù)性，采用“干濕水深處理方式”，將干、濕水深及淹沒水深分別設置為0.005m、0.01m與0.05m，渦黏系數(shù)根據(jù)Smagorinsky公式設置為0.28。

2.3 一、二維水動力模型的耦合

本文采用MIKE FLOOD將MIKE11和MIKE21 FM進行耦合，耦合時采用側向連接來模擬洪水從河道漫流至洪泛區(qū)的運動過程。連接方式允許MIKE21 FM的一系列網(wǎng)格單元從側面連接到MIKE11的部分河段甚至整個河段[11]。連接方法采用“Cell to Cell”，即對結構物中每1個節(jié)點都進行水動量方程計算，之后所有經(jīng)過計算的水流都被重新分配到MIKE11和MIKE21 FM網(wǎng)格單元中[15]。

3 模型糙率的率定與驗證

3.1 率定

3.1.1一維模型糙率的率定

選取1998年6—10月汛期實測水位及沿程洪痕數(shù)據(jù)用于河道糙率的率定，根據(jù)文獻[16]及靳宏偉等[17-18]的松花江流域中小洪水糙率，并結合河道特性與現(xiàn)場查勘情況，確定模型參數(shù)，得出河道下游出流過程，然后提取河段沿程各斷面水位峰值與實測洪痕值進行對比，并參考中間站白沙灘的水位擬合情況進行綜合分析。經(jīng)率定，確定模型計算步長為5 s，河段綜合糙率為0.042 7，該糙率符合《松花江流域防洪規(guī)劃》。

模型模擬結果與實測結果對比如圖2所示。由圖2(a)可知，嫩江江橋—大賚河段沿程各斷面的水位峰值誤差均小于0.2 m，該結果符合GB/T22482—2008《水文情報預報規(guī)范》[19]評定要求。中間站白沙灘的納什效率系數(shù)(NSE)[20]、決定系數(shù)(R2)、一致性指數(shù)(d)與均方根誤差(RMSE)[21-23]分別為0.98、0.98、0.99和0.174，說明模型擬合結果較好，模擬水位變化過程與實測水位變化過程基本一致，如圖2(b)所示。

圖2 1998年江橋—大賚河段洪痕率定

3.1.2一、二維耦合模型糙率的率定

由于嫩江流域洪水實際淹沒范圍矢量圖難以獲取，本文基于USGS官網(wǎng)(https://earthexplorer.usgs.gov/)下載了洪水時期研究區(qū)的Landsat影像，包括Landsat 7和Landsat 8。然后利用ENVI遙感圖像處理方法中的監(jiān)督分類選取洪水淹沒區(qū)、未淹沒區(qū)與其他樣本，并通過支持向量機分類器進行分類，然后結合遙感目視解譯提取的洪水實際淹沒范圍與模擬的淹沒范圍進行對比?？紤]到汛期遙感影像圖受云層影響較大，本文選取2013年8月11日(Landsat 8,云量低于20%)作為一、二維耦合模型糙率的率定時段。率定過程中，模型上邊界選用江橋站2013年8月11日前后的實測流量過程，下邊界選取大賚站同期的水位過程。經(jīng)率定，嫩江流域耕地、林地、草地、濕地、水體、人工表面及裸地糙率分別為0.06、0.065、0.04、0.035、0.021、0.07和0.035，時間步長為30 s。模擬淹沒范圍與實際淹沒范圍對比結果如圖3所示。由圖3(a)可知，2013年8月11日洪水實際淹沒范圍為764.50 km2，模型模擬淹沒范圍為761.72 km2，模擬與實際相交的淹沒范圍為700.83 km2，模型相對誤差為0.4%，滿足《洪水風險圖編制導則》[24]的要求。

圖3 嫩江江橋—大賚河段率定與驗證結果

3.2 糙率的驗證

3.2.1一維模型糙率的驗證

為保證江橋—大賚河段一維水動力模型的計算精度及普適性，選取2013年典型洪水年的實測水文數(shù)據(jù)進行驗證，驗證結果如圖4所示。

圖4 2013年白沙灘水位模擬與實測結果對比

由圖4可知，2013年白沙灘實測水位峰值為136.49 m，模擬水位峰值為136.483 m，絕對誤差為0.007 m，相對誤差為0.01%，NSE、R2、d與RMSE分別為0.97、0.98、0.99和0.192。由此可見，2013年實測水位峰值與模擬水位峰值的絕對誤差與相對誤差均滿足《水文情報預報規(guī)范》的峰值誤差精度要求；納什效率系數(shù)、測定系數(shù)、一致性指數(shù)與均方根誤差離散程度低、擬合效果好，即模擬水位變化過程與實測過程基本一致，其中納什效率系數(shù)符合《水文情報預報規(guī)范》中的甲級預報精度標準。因此，該模型選用的綜合糙率模擬精度高，可為嫩江流域的洪水預報、防洪決策及生態(tài)河道建設提供技術支撐。

3.2.2一、二維耦合模型糙率的驗證

考慮到汛期遙感影像圖受云層影響較大，選取2013年9月4日(Landsat 7，云量低于20%)作為一、二維耦合模型糙率的驗證時段。驗證過程中，模型上邊界選用江橋站2013年9月4日前后的實測流量過程，下邊界選取大賚站同期的水位過程，模型模擬的淹沒范圍與實際淹沒范圍對比如圖3(b)所示。由圖3(b)可知，2013年9月4日提取的實際淹沒范圍為942.70 km2，模擬淹沒范圍為934.49 km2，模擬與實際相交的淹沒范圍為895.83 km2，模型相對誤差為 1%，滿足《洪水風險圖編制導則》的要求。

綜上所述，本文在一維模型和一、二維耦合模型計算中，選用江橋—大賚河段多個實測斷面的洪痕值進行模型糙率的率定與驗證，且率定與驗證的持續(xù)時間較長(時間范圍為6—10月)，這樣可有效避免選取單一水文站進行糙率率定與驗證的偶然性，同時也采用了中間站白沙灘的實測水位過程進行擬合對比，有效增加了糙率率定與驗證的合理性。

4 不同重現(xiàn)期洪水演進及其生態(tài)學意義

4.1 不同重現(xiàn)期洪水淹沒范圍數(shù)值模擬

為了更好地促進嫩江流域生態(tài)河道建設及防洪決策的制定，本文就不同重現(xiàn)期洪水淹沒范圍進行模擬。其中上邊界的流量過程主要參考《松花江流域防洪規(guī)劃》中的洪峰流量進行放大，下邊界的控制條件選取大賚站的水位-流量關系。利用MIKE FLOOD一、二維水動力耦合模型進行不同頻率洪水淹沒范圍數(shù)值模擬，結果如圖5所示。圖5(a)為不同重現(xiàn)期洪水淹沒范圍，圖5(b)為典型斷面示意圖(斷面起點位置為124°3.740′E、45°43.906′N，終點位置為124°8.349′E、45°47.197′N)。

圖5 嫩江江橋—大賚河段不同重現(xiàn)期洪水淹沒范圍

從模型計算水位分析，在5年一遇洪水模擬工況下，洪水演進至8月6日后(白沙灘附近洪水由主河槽溢出灘地)，江橋—大賚河段總淹沒面積為58.50 km2，最大水位高程為138.09 m。洪水演進至8月22日后，江橋到大賚河段淹沒范圍最大，此時最大淹沒范圍為677.14 km2，最大水深為12.09 m。

在20年一遇洪水模擬工況下，洪水演進至8月6日后，江橋—大賚河段總淹沒面積為192.69 km2,最大水位高程為139.75 m。洪水演進至8月22日后，江橋到大賚河段淹沒范圍最大，此時最大淹沒范圍為896.06 km2，最大水深為13.86 m。

在50年一遇洪水模擬工況下，洪水演進至8月6日后，江橋—大賚河段總淹沒面積為230.46 km2,最大水位高程為140.48 m。洪水演進至8月22日后，江橋到大賚河段淹沒范圍最大，此時最大淹沒范圍為1 044.26 km2，最大水深為15.05 m。

在100年一遇洪水模擬工況下，洪水演進至8月6日后，江橋—大賚河段總淹沒面積為261.58 km2,最大水位高程為141.16 m。洪水演進至8月22日后，江橋到大賚河段淹沒范圍最大，此時最大淹沒范圍為1 234.43 km2，最大水深為15.53 m。

由此可知，隨著洪水頻率的增加，嫩江江橋—大賚河段的洪水淹沒范圍也不斷增加(圖5)，其中20年一遇較5年一遇洪泛增加面積為218.92 km2，50年一遇較20年一遇洪泛增加面積為148.20 km2，100年一遇較50年一遇洪泛增加面積為190.17 km2。

4.2 不同重現(xiàn)期洪水擾動對洪泛區(qū)植被指數(shù)的影響

洪泛區(qū)位于河流與陸地環(huán)境的過渡地帶，其植被的生長需要一定范圍的洪水淹沒來提供充足的水分[7]，但洪水淹沒的頻繁程度又將在一定程度上影響植被的正常生長。由于植被具有不同的光譜特征規(guī)律，目前在遙感應用領域?qū)⒖梢姽馀c近紅外遙感光譜觀測通道進行組合運算, 便可得到衡量地表植被狀況的簡單、有效和經(jīng)驗度量的植被指數(shù)，其中歸一化植被指數(shù)[25](normalized difference vegetation index，NDVI)是由紅光、近紅外2個波段計算生成的，其基本公式為

(1)

式中：ρNIR為近紅外波段的表觀反射率或方向放射率；ρRED為紅光波段的表觀反射率或方向放射率。

增強型植被指數(shù)(enhanced vegetation index，EVI) 是由Red、NIR、Blue(紅光、近紅外、藍光)3個波段計算生成的，其基本公式為

(2)

式中:ρBLUE為藍光波段的表觀反射率或者方向放射率；G為增益系數(shù);C1、C2為大氣調(diào)節(jié)參數(shù)；L為土壤調(diào)節(jié)參數(shù)。在EVI計算中，通常取L=1，C1=6，C2=7.5，G=2.5。上述兩種植被指數(shù)作為常用的植被指數(shù)，已廣泛用來定性和定量評價植被覆蓋及其生長活力[26]，對于指示植被的生長狀態(tài)與空間分布具有重要意義。

4.2.1洪泛區(qū)植被指數(shù)季節(jié)變化分析

為深入了解嫩江流域植被指數(shù)的季節(jié)變化規(guī)律，本文利用ArcGIS對不同洪水頻率淹沒范圍內(nèi)2000—2019年的MOD13Q1 植被指數(shù)數(shù)據(jù)(時間分辨率16 d，空間分辨率250 m，全年共23景)進行掩膜[7]，并對每一景影像求像元平均值，分析結果如圖6所示。由圖6可見，2000—2020年不同重現(xiàn)期內(nèi)的植被指數(shù)(NDVI、EVI)具有較強的季節(jié)性變化規(guī)律，6—8月夏季植被生長最為旺盛，7月底至8月初植被指數(shù)達到最大值，隨后受洪水影響植被指數(shù)從峰值急劇下降，8—9月降低幅度最大，同時受物候條件影響植被在10月底降至完全枯黃，植被指數(shù)進入“低谷階段”，且不同洪水頻率下的NDVI和EVI的變化一致性較好。該結論與饒品增等[6]的2013年黑龍江省洪水對植被的影響評估結果一致。

圖6 2000—2020年植被指數(shù)變化

圖7 2012—2016年植被指數(shù)變化

4.2.2典型洪水年前后植被指數(shù)變化分析

選取2013典型洪水年前后的植被指數(shù)進行分析(圖7)，該年份洪水重現(xiàn)期達10年一遇。由圖7可知，在2013典型洪水年前后，不同重現(xiàn)期內(nèi)的植被指數(shù)(NDVI、EVI)存在顯著差異，2013年與2012年相比，2013年不同重現(xiàn)期內(nèi)的NDVI、EVI峰值均明顯降低，可見洪水對該區(qū)域的植被造成了顯著的影響，其中5年一遇植被指數(shù)降低幅度最大，20年一遇較5年一遇洪泛增加區(qū)的植被指數(shù)也呈現(xiàn)一定的下降趨勢；與2013年對比，2014年不同重現(xiàn)期內(nèi)的NDVI和EVI的峰值均有較明顯的提高，其中5年一遇的植被指數(shù)增加最多，但與2012年的植被指數(shù)峰值還存在一定差距，可見，2014年由于洪水干擾強度降低，洪泛區(qū)的植被指數(shù)均得到了明顯的回升，但回升幅度尚未完全達到洪水前狀況，而2015年洪泛區(qū)內(nèi)的植被指數(shù)回升則逐漸接近洪水前狀況。由此可見，不同重現(xiàn)期洪水對植被具有明顯的破壞作用，但這種破壞作用在嫩江流域周期性相對較短，表明該地區(qū)植被的可恢復性較強，植被類型主要以1～2 a生為主，與饒品增等[6]的研究結果一致。

4.2.3不同重現(xiàn)期洪泛區(qū)植被指數(shù)變化分析

為探究不同重現(xiàn)期洪泛區(qū)內(nèi)植被指數(shù)的變化規(guī)律，本文統(tǒng)計了2000—2020年植被指數(shù)的年平均變化值，如圖8所示。不同重現(xiàn)期洪泛區(qū)內(nèi)的植被指數(shù)存在顯著差距，其中頻繁淹沒區(qū)(5年一遇以下)的植被指數(shù)最小，中等頻繁淹沒區(qū)(50年一遇較5年一遇洪泛增加區(qū))的植被指數(shù)顯著增大，50年一遇較20年一遇洪泛增加區(qū)的植被指數(shù)增加達到最大，而較低頻繁淹沒區(qū)(100年一遇較50年一遇洪泛增加區(qū))的植被指數(shù)出現(xiàn)下降趨勢。這說明不同重現(xiàn)期洪水擾動對洪泛區(qū)植被指數(shù)產(chǎn)生了一定程度的影響，該結果與加拿大魁北克南部不同洪水頻率對物種多樣性的影響結果趨于一致，同時也符合Connell[27]的中度干擾假說。

圖8 2000—2020年植被指數(shù)的年際變化

根據(jù)上述洪水演進數(shù)值模擬及植被指數(shù)數(shù)據(jù)分析發(fā)現(xiàn)，利用不同重現(xiàn)期的洪水淹沒范圍可初步對嫩江流域洪泛區(qū)植被群落進行分區(qū)，這不僅為相關的野外植被調(diào)查提供了前期準備，也為嫩江流域防洪規(guī)劃的制定及生態(tài)河道建設提供了基礎數(shù)據(jù)與技術支撐。

5 結 論

a.嫩江下游一、二維水動力耦合模型在率定期與驗證期的模擬效果符合甲級預報精度，模型模擬精度較高。

b.不同重現(xiàn)期內(nèi)的植被指數(shù)具有較強的季節(jié)變化規(guī)律，夏季植被生長的最為旺盛，且NDVI與EVI的變化基本一致。同時洪水對植被具有明顯的破壞作用，但這種破壞作用在嫩江流域周期性相對較短，表明該地區(qū)植被的可恢復性較強，植被類型主要以1～2 a生為主。

c.受洪水頻繁淹沒程度影響，頻繁淹沒區(qū)的植被指數(shù)最小，中等頻繁淹沒區(qū)的植被指數(shù)顯著增大，50年一遇較20年一遇洪泛增加區(qū)內(nèi)的植被指數(shù)達到最大，而較低頻繁淹沒區(qū)的植被指數(shù)則呈現(xiàn)下降趨勢，該結論符合中度干擾假說。