關(guān)銅壘,劉佳嘉,周祖昊,井 涌,王鵬翔,杜 崇
(1.黑龍江大學 水利電力學院,黑龍江 哈爾濱 150080;2.中國水利水電科學研究院流域水循環(huán)模擬與調(diào)控重點實驗室,北京 100038;3.陜西省水文水資源勘測局,陜西 西安 710068)
分布式水文模型計算單元劃分方法主要有網(wǎng)格單元、地貌單元和子流域單元[1]3 種,其中:基于坡面漫流模型的子流域單元劃分方法應用較為廣泛且被成功集成到各類GIS 軟件及水文分析軟件(如ArcGIS、ARCSWAT、QSWAT 等)中[2-3]。一些學者在前人研究基礎上進行改進,以滿足不同情形的子流域劃分要求,如劉歡等[2]基于傳統(tǒng)子流域劃分方法,提出一種面向大尺度區(qū)域的改進方法,可實現(xiàn)流域出水口的自動識別,并通過多閾值模擬河網(wǎng)融合技術(shù),解決河網(wǎng)提取失真問題所導致的子流域拓撲關(guān)系錯亂問題;雷曉輝等[4-5]提出一種子流域劃分過程中界河、海岸線處理方法,并基于傳統(tǒng)子流域劃分方法提出一種改進的子流域加密算法,可解決狹長子流域和水庫、水文站控制的子流域劃分問題;劉佳嘉等[6]提出一種能夠描述城市管網(wǎng)匯流作用的子流域劃分方法,可有效解決城市管網(wǎng)對子流域劃分的影響。然而無論是傳統(tǒng)子流域劃分方法還是改進后的方法,在進行子流域劃分時都未考慮流域內(nèi)湖泊的相對獨立性,一般通過虛擬河網(wǎng)將湖泊劃分為若干部分,每個部分分屬不同的子流域。湖泊作為流域的重要組成部分,其水量平衡過程對流域水循環(huán)過程的影響不可忽視,現(xiàn)有劃分方法使得湖泊水量平衡計算變得比較困難,因此有必要將湖泊劃分在同一個子流域內(nèi)進行研究分析,但在將湖泊進行獨立劃分的時候,會遇到湖泊上游河流匯入和坡面區(qū)域編碼等問題。鑒于此,筆者提出一種考慮湖泊范圍的子流域劃分和編碼方法,對湖泊整體及其上游匯入?yún)^(qū)域進行準確劃分及編碼,并利用流域分布式水文模型(WEP-L 模型)[7-8]進行模擬驗證。
選取黃河源區(qū)(唐乃亥以上區(qū)域)作為研究區(qū)。研究區(qū)總面積12.2 萬km2,包含黃河源頭兩大淡水湖泊扎陵湖和鄂陵湖。這兩個湖泊均位于瑪多縣,其中:扎陵湖東西長35 km、南北寬21.6 km,湖面面積526 km2,湖泊東北部較深、西部較淺,平均水深約9 m,蓄水量46.7 億m3;鄂陵湖位于扎陵湖之東,南北長32.3 km、東西寬31.6 km,湖面面積610 km2,平均水深約18 m,蓄水量107 億m3[9]。
本研究所用DEM(見圖1)來源于中國科學院信息中心提供的ASTER GDEMV2 高分辨率數(shù)字高程數(shù)據(jù),空間分辨率為30 m×30 m。為提高計算效率,采用ArcGIS 軟件將該數(shù)據(jù)重采樣為1 km×1 km 分辨率數(shù)據(jù)。流域水系數(shù)據(jù)來源于國家基礎地理信息中心的全國1 ∶250 000 地形數(shù)據(jù)庫。
圖1 研究區(qū)DEM
模擬河網(wǎng)提取過程采用傳統(tǒng)方法進行,即經(jīng)填洼、流向計算、匯流累計數(shù)計算、河道閾值設置、模擬河網(wǎng)提取等操作得到[10-11]。首先,根據(jù)河網(wǎng)和湖泊對流域柵格屬性進行設置,主要包括A和B兩部分,其中:A=1 代表湖泊、A=0 代表非湖泊;B=1 代表河道、B=0 代表非河道(坡面)。因此,流域柵格屬性有湖泊河道(AB=11)、湖泊非河道(AB=10)、非湖泊河道(AB=01)以及非湖泊非河道(AB=00)4 種。
根據(jù)湖泊范圍及匯流關(guān)系,將流域分為湖泊區(qū)域、湖濱帶區(qū)域、其他區(qū)域3 類區(qū)域,其中:湖泊區(qū)域為湖泊范圍內(nèi)部區(qū)域,湖濱帶區(qū)域為直接匯入湖泊的無河網(wǎng)坡面集水區(qū)域,其他區(qū)域則是除湖泊和湖濱帶以外的區(qū)域。湖泊區(qū)域和湖濱帶區(qū)域采用湖區(qū)子流域劃分方法進行劃分,其他區(qū)域采用非湖區(qū)子流域劃分方法進行劃分。兩種方法的選擇主要由流域柵格屬性確定,即在對河段從下游往上游追溯劃分的時候,如果該河段起始河道柵格為非湖泊河道柵格(AB=01),則采用非湖區(qū)子流域劃分方法進行劃分;如果該河段起始河道柵格為湖泊河道柵格(AB=11),則采用湖區(qū)子流域劃分方法進行劃分。
3.2.1 非湖區(qū)子流域劃分方法
非湖區(qū)子流域劃分方法即傳統(tǒng)子流域劃分方法,具體步驟[5,11-13]:①從流域出口柵格開始,基于柵格流向沿河道從下游向上游追溯,分河段進行編碼;②當完成一個子流域劃分后,沿河道向上對該河段的上游河道繼續(xù)追溯,根據(jù)流入河道柵格的匯流量確定干支流,將匯流量最大的河道作為干流,將其余的河道作為支流;③先對匯流量小的支流河道進行子流域劃分,再對匯流量大的干流河道進行子流域劃分,直至河流源頭或者遇到湖泊河道柵格結(jié)束。
3.2.2 湖區(qū)子流域劃分方法
湖區(qū)子流域劃分方法是本研究的創(chuàng)新點,主要思路是將每個湖泊劃分為一個獨立的湖泊單元,并且根據(jù)河網(wǎng)與湖泊進出口節(jié)點將匯入該湖泊的湖濱帶區(qū)域劃分為若干個單元。具體步驟如下。
(1)湖區(qū)區(qū)域預處理。以湖泊范圍內(nèi)的模擬河網(wǎng)與湖泊出入口節(jié)點為分割點,對湖泊進行分割處理,將兩條模擬河網(wǎng)之間的湖泊范圍內(nèi)柵格作為一個整體,繼而將湖泊范圍分割為若干個子區(qū)域,每一個區(qū)域內(nèi)的柵格相連通且由模擬河網(wǎng)分割開。對分割開的子區(qū)域賦予不同的分區(qū)編號X,該編號是從1 開始的連續(xù)自然數(shù)(如圖2 的1~4 號),流域內(nèi)的每個湖泊分割后的子區(qū)域編號均從1 開始。
圖2 湖泊預處理結(jié)果示意
(2)湖泊單元劃分。當湖泊下游子流域劃分編碼完成后,溯源到當前湖泊出口柵格時開始湖區(qū)子流域劃分。從該柵格出口開始溯源編碼,直至遇到非湖泊柵格結(jié)束(AB=00 或AB=01),將其中所有湖泊柵格(AB=10 或AB=11)都賦予相同的編號,如圖3 的6 號和12 號湖泊型子流域。
圖3 子流域劃分結(jié)果示意
(3)湖濱帶單元劃分。根據(jù)湖泊預處理結(jié)果,以湖泊范圍內(nèi)擁有不同分區(qū)編號的非河網(wǎng)柵格為起點,對湖濱帶柵格進行溯源編號,將湖濱帶流入湖泊內(nèi)不同分區(qū)的湖濱帶柵格賦予相應的子流域編號N+x(N為湖濱帶柵格流入的湖泊單元編號,x 為湖濱帶柵格最終流入的湖泊柵格所擁有的湖泊子區(qū)域編號,如圖3 中7~9 號子流域分別對應6+1、6+2、6+3)。重復上述步驟,直至該湖泊上游所有湖濱帶柵格都完成子流域劃分,至此完成當前的湖區(qū)子流域劃分,然后根據(jù)湖泊上游匯入河道柵格累計數(shù)按從小到大的順序依次執(zhí)行非湖區(qū)子流域劃分。
3.2.3 確定子流域類別屬性
所有子流域劃分完成后,根據(jù)子流域內(nèi)的湖泊范圍和河道情況對子流域進行屬性劃分,主要包括一般類型、湖泊類型和坡面類型3 種。其中:一般類型的子流域包括一條河道柵格以及連接該河道的坡面柵格,無湖泊柵格(如圖3 的1~2 號子流域);湖泊類型的子流域包含該湖泊范圍內(nèi)的所有河道和坡面柵格(如圖3 的6 號、12 號子流域);坡面類型的子流域為連接湖泊區(qū)域的坡面柵格(如圖3 的7~8 號、13~14 號子流域),不含河道和湖泊柵格。
3.3.1 確定上游流入子流域
本研究中由于存在不包括河網(wǎng)的坡面類型的子流域,無法直接根據(jù)河網(wǎng)匯流次序確定子流域間上下游拓撲關(guān)系,因此需要根據(jù)柵格流向進行溯源遍歷(沿模擬河網(wǎng)對流域內(nèi)所有柵格逐個訪問)確定,具體規(guī)則及步驟如下。
(1)如果遍歷的子流域為一般類型子流域,則根據(jù)子流域內(nèi)的河道柵格進行溯源遍歷,查找上游流入的河道所在子流域作為其上游(如圖3 一般類型1 號子流域的上游2 號和5 號子流域、一般類型5 號子流域的上游6 號子流域)。
(2)如果遍歷的子流域為湖泊類型子流域,則需要先對上游匯入的湖濱帶子流域進行遍歷,再對匯入的河道所在子流域進行遍歷。首先,根據(jù)柵格流向,查找湖濱帶直接流入湖泊的坡面類型子流域,根據(jù)湖泊預處理結(jié)果,將湖濱帶直接流入湖泊不同子區(qū)域的多個湖濱帶子流域都作為湖泊類型子流域的上游坡面類型子流域(如圖3 湖泊類型6 號子流域的上游7~9 號子流域、12 號子流域的上游13~16 號子流域);其次,根據(jù)河道柵格流向,沿湖泊內(nèi)虛擬河道向上繼續(xù)追溯,將匯入的非湖泊河道柵格所在的子流域作為當前湖泊類型子流域的上游(如圖3 湖泊類型6 號子流域的上游11號子流域、12 號子流域的上游17~19 號子流域)。
(3)如果遍歷的子流域為坡面類型子流域,則直接結(jié)束遍歷,無上游子流域,如圖3 中坡面類型7 號子流域。
3.3.2 建立上下游拓撲關(guān)系屬性表
在確定子流域上游流入子流域后,需要構(gòu)建子流域上下游拓撲關(guān)系屬性表,上下游拓撲關(guān)系屬性表內(nèi)容主要包括當前子流域編號、下游子流域編號、上游子流域數(shù)目、多個上游子流域編號。
按照上述方法對黃河源區(qū)進行子流域劃分。首先采用ArcGIS 軟件經(jīng)過柵格填洼、計算柵格流向、計算柵格匯流量、確定河網(wǎng)閾值等步驟提取模擬河網(wǎng),然后基于提取的模擬河網(wǎng)設置流域柵格屬性,接著根據(jù)鄂陵湖和扎陵湖邊界及模擬河網(wǎng)對兩個湖泊進行預處理,最后對研究區(qū)域進行子流域劃分以及構(gòu)建上下游拓撲關(guān)系屬性表。黃河源區(qū)共劃分1 871 個子流域,將兩個湖泊劃分為單獨的子流域,其中:鄂陵湖上游有河道匯入的子流域4 個、無河道匯入的湖濱帶坡面型子流域5 個,扎陵湖有河道匯入的子流域5 個、無河道匯入的湖濱帶坡面型子流域7 個。計算單元劃分結(jié)果如圖4 所示,兩個湖區(qū)子流域的上下游拓撲關(guān)系屬性表見表1。
表1 湖區(qū)子流域上下游拓撲關(guān)系屬性
圖4 黃河源區(qū)考慮湖泊范圍的子流域劃分結(jié)果
為對比本文提出的考慮湖泊范圍的子流域劃分方法與傳統(tǒng)子流域劃分方法之間的優(yōu)缺點,采用ArcGIS軟件對黃河源區(qū)(唐乃亥以上區(qū)域)進行子流域劃分,結(jié)果見圖5(左下角為湖泊范圍子流域劃分放大圖)。對比圖4 和圖5 可以發(fā)現(xiàn),兩種方法非湖泊范圍劃分的子流域大體相當,但傳統(tǒng)方法將湖泊范圍劃分為若干個子流域,湖泊范圍中心處出現(xiàn)密集的異常子流域;而本文提出的方法將湖泊范圍單獨劃分為一個子流域,保證了湖泊的獨立性。
圖5 黃河源區(qū)常規(guī)方法子流域劃分結(jié)果
采用流域分布式水文模型WEP-L 進行模擬驗證。為符合本文提出的子流域劃分單元的特點,對該模型的匯流模塊加以改進,即在匯流模塊中識別湖濱帶計算單元,因其內(nèi)部不包含河網(wǎng),故只對其進行坡面匯流演算,而不進行河道匯流計算。選取瑪曲和唐乃亥水文站1956—2000 年月均流量進行水文模擬驗證,其中1956—1980 年為率定期、1981—2000 年為驗證期,并將Nash 效率系數(shù)和相對誤差作為水文模擬結(jié)果的檢驗指標,結(jié)果如圖6 所示。瑪曲水文站1956—1980 年率定期的Nash 效率系數(shù)為0.764、相對誤差為1.6%,1981—2000 年驗證期的Nash 效率系數(shù)為0.750、相對誤差為-2.3%;唐乃亥水文站1956—1980 年率定期的Nash 效率系數(shù)為0.806、相對誤差為2.1%,1981—2000 年驗證期的Nash 效率系數(shù)為0.768、相對誤差為-1.4%。月徑流量模擬結(jié)果較好,表明本文提出的子流域劃分方法符合分布式水文模型建模的需求,能有效模擬水文站月徑流量過程。本文所提出的方法優(yōu)勢在于能通過水文模擬方便地平衡湖泊的水量,為進一步研究湖泊水循環(huán)及伴生生態(tài)演變提供支持。在改進計算單元劃分的前提下,借助WEP-L 模型輸出鄂陵湖和扎陵湖的蒸散發(fā)量(見圖7)和出、入湖水量(見圖8)。扎陵湖的年均蒸散發(fā)量為597.7 mm,鄂陵湖的年均蒸散發(fā)量為746.9 mm;鄂陵湖多年平均出、入湖水量分別為16.1 億、17.2 億m3,扎陵湖多年平均出、入湖水量分別為8.3 億、9.8 億m3。
圖6 瑪曲和唐乃亥水文站月徑流量模擬結(jié)果
圖7 鄂陵湖和扎陵湖多年平均蒸散發(fā)量
圖8 鄂陵湖和扎陵湖多年平均出入湖水量
(1)本文提出一種考慮湖泊范圍的子流域劃分和編碼方法,該方法通過對湖泊范圍進行預處理,可自動區(qū)分湖泊類型子流域上游匯入的河道和坡面區(qū)域,將上游匯入的坡面湖濱帶劃分為若干個子流域,實現(xiàn)對研究區(qū)內(nèi)湖泊區(qū)域的自動識別,并將湖泊劃分為單獨的子流域然后編碼。
(2)與現(xiàn)有子流域劃分方法相比,本文提出的改進方法可將湖泊劃分為單獨的子流域,保證湖泊的完整性,便于運用分布式水文模型計算湖泊的進、出水量。
(3)采用流域分布式水文模型WEP-L 進行月均流量模擬驗證,結(jié)果表明所提出的子流域劃分方法能夠滿足分布式水文模型建模的需求,保證子流域劃分過程中湖泊的獨立性、完整性和有效性。