基于EnKS和SWAT模型的閩江流域徑流數據同化

項 勇, 陳蕓芝, 唐麗芳, 汪小欽

(1.福州大學 數字中國研究院(福建), 福建 福州 350108; 2.福州大學 空間數據挖掘和信息共享教育部重點實驗室,福建 福州 350108; 3.福建省水土保持實驗站, 福建 福州 350001)

1 研究背景

流域水文模型可以反映整個流域的綜合水文信息,在水旱災害防治、水資源管理和開發利用等方面應用廣泛[1-3]。但受模型結構、觀測數據、氣象輸入等不確定性因素的影響,水文模型在進行水文狀態模擬時存在不可避免的偏差[4]。數據同化技術能夠綜合不同觀測信息,以獲得流域狀態變量的最佳估計,從而減少水文模型徑流模擬偏差[5]。已有研究表明,同化遙感數據,如土壤濕度、蒸散發、積雪等可以有效減少模型輸出的不確定性[6-8]。然而由于徑流觀測數據與徑流模擬數據直接相關,且地表實測徑流數據具有較高可靠性,同化站點徑流實測信息仍是徑流預測的最優選擇[9]。

數據同化方法可分為變分數據同化和順序數據同化。集合卡爾曼濾波(ensemble Kalman filter, EnKF)是一種順序同化方法,可以根據實測信息順序更新模型狀態。相較于普通卡爾曼濾波,EnKF解決了實際應用中背景誤差協方差矩陣計算困難的問題,并且能夠直接應用于非線性模型[10]。同時,相較于更復雜的粒子濾波等方法,EnKF的計算效率更高,在水文數據同化中得到了廣泛的應用[11-13]。標準EnKF方法在小流域尺度以及狀態變量與觀測變量之間實時響應的情況下取得了良好的效果[14-15],但該方法未考慮土壤水等水文狀態變量和河道出口流量之間的時間滯后問題,導致中大型流域的徑流數據同化效果較差[16]。針對時間滯后問題,學者們提出了不同解決方法,如Weerts等[17]設置了固定滯后時間,通過當前觀測值更新固定滯后時間前的模型狀態變量;McMillan等[18]基于REnKF方法,利用當前觀測值遞歸地更新最大滯后期前至當前觀測時間的模型狀態變量,結果表明該方法可以解決集水區土壤含水量和觀測流量之間的響應滯后問題,但其使用同一觀測數據迭代更新先前模型狀態,存在過度校正的風險,且計算復雜。集合卡爾曼平滑器(ensemble Kalman smoother, EnKS)是EnKF的重要變體,相較于EnKF方法,EnKS綜合考慮了時間窗口內的狀態變量和觀測值,使得預測結果更為平滑、準確,能夠有效解決徑流數據同化中的時間滯后問題。該方法在集總式水文模型中表現良好[10,19],但目前仍較少應用于分布式和半分布式水文模型。

本文以閩江流域為研究區,基于EnKS方法和半分布式水文模型SWAT(soil and water assessment tool),構建徑流數據同化方案,在確定最佳時間窗口長度的前提下,評估EnKS方法對半分布式水文模型的優化效果,分析數據同化對于不同徑流分量的影響,為閩江流域徑流數據同化提供參考。

2 數據來源與研究方法

2.1 研究區概況

本文選取閩江流域為研究區,如圖1所示。研究區流域面積為59 922 km2,地處東經116°23′～119°35′,北緯25°23′～28°16′,主要支流包括建溪、富屯溪和沙溪等。閩江流域屬于亞熱帶季風氣候,雨量充沛,但分布不均,降水主要集中在夏秋季節,春冬季節降水較少,多年平均降水量為1 710 mm。

圖1 閩江流域水系及水文、氣象、雨量站點分布

2.2 數據來源

閩江流域SWAT模型構建所使用的數據如表1所示。氣象數據為中國氣象數據網(http://data.cma.cn/)提供的泰寧、福州、浦城、建甌等8個氣象站點2017—2020年逐日氣溫、相對濕度、風速等數據,以及福建省水文水資源勘測局提供的寧化、大田、松溪、武夷山等10個雨量站點2017—2020年逐日降水數據。實測徑流數據為福建省水文水資源勘測局提供的七里街水文站、沙縣水文站、竹岐水文站2018—2020年逐日實測徑流量數據。

表1 閩江流域SWAT模型構建中使用的數據

2.3 SWAT模型

SWAT模型是應用最廣泛的半分布式水文模型之一[20],該模型根據地形將整個流域劃分為多個子流域,并根據下墊面條件進一步劃分水文響應單元(hydrological response unit,HRU),HRU是SWAT最基本的計算單元。

SWAT模擬的逐日水文循環基于每個HRU的水量平衡方程進行,其計算方程如下:

(1)

式中:SWt為第t天結束時的土壤含水量,mm;SW0為初始土壤含水量,mm;Ri為第i天的總降水量,mm;Qsurf,i為第i天的地表徑流量,mm;Ei為第i天的蒸散量,mm;Wseep,i為第i天的下滲量和壤中流,mm;Qgw,i為第i天的地下水徑流量,mm。

總徑流Qi由第i天的地表徑流Qsurf,i、壤中流Qlat,i和地下水徑流Qgw,i組成:

Qi=Qsurf,i+Qlat,i+Qgw,i

(2)

各HRU產生的總徑流在子流域尺度上匯集,繼而進入河網匯總,并在河道出口斷面處計算徑流。

2.4 集合卡爾曼濾波

集合卡爾曼濾波基于Monte Carlo方法和標準卡爾曼濾波公式來估計非線性模型中的狀態和參數[10]。EnKF主要分為預報和更新兩步。

在預報步驟中,將初始化后的k時刻狀態變量集合代入模型,得到k+1時刻的狀態變量預測值:

(3)

在基于模型狀態預報值獲得模型狀態預報協方差矩陣后,模型狀態變量的預報值基于觀測信息更新,其公式如下:

(4)

(5)

(6)

(vk+1～N(0,Rk+1))

(7)

2.5 集合卡爾曼平滑器

EnKS是EnKF的擴展,該方法將平滑窗口內存在的全部觀測與模擬分別納入觀測變量值集合與模型狀態變量集合,然后通過EnKF公式更新狀態變量。EnKF方法即為EnKS時間窗口等于1天的特殊情況。

在初始化隨機變量集合后,SWAT模擬運行至平滑窗口最終時刻,獲取窗口內所有狀態變量模擬值,因此擴展的狀態變量X包括平滑窗口內所有L個狀態值x,其可表示為:

X=[x1x2…xL]T

(8)

相應地,擴展的觀測向量集合Y包括時空獨立分布的所有m個觀測值,即為:

Y=[y1y2…ym]T

(9)

而后使用EnKF公式計算卡爾曼增益矩陣,對模型狀態變量進行更新。

2.6 評價指標

為評估模型預測與實測值的接近程度,采用納什效率系數(Nash-Sutcliffe efficiency coefficient,NSE)、均方根誤差(root mean square error,RMSE)共同評估模型模擬效果[21]。其中,NSE反映了模擬值和實測值的擬合度,NSE越趨近于1則模擬效果越好;RMSE反映了模擬值與實測值的偏差,RMSE越趨近于0,代表模擬效果越好。

3 徑流數據同化實驗

圖2 徑流觀測數據同化方案流程

3.1 SWAT模型參數率定

選取閩江流域七里街、沙縣、竹岐3個水文站實測徑流量數據對SWAT模型進行率定和驗證,以2018年1月1日—2019年12月31日為率定期,2020年1月1日—2020年12月31日為驗證期,流域參數率定結果如表2所示。

表2 流域參數率定結果

結果表明,七里街、沙縣、竹岐3個水文站的NSE在率定期分別為0.76、0.74、0.82,驗證期分別為0.70、0.50、0.72。一般認為當NSE>0.5時,模型模擬結果可被接受[22]。因此該SWAT模型可用于后續閩江流域徑流數據同化研究。

3.2 狀態變量選擇與更新

土壤水決定地表水文循環的諸多過程,是水文模型模擬中最重要的變量之一,在SWAT同化中可以將所有HRU的土壤含水量納入狀態向量中進行計算,但這增加了同化的計算量,因此本文定義了土壤水校正比率分別控制各個站點上游的土壤含水量,共設置3個比率,其控制范圍如圖3所示。HRU尺度上各層土壤含水量更新公式如下:

圖3 研究區不同土壤水校正比率空間分布

sol_sta(c,b)=sol_st(c,b)·ratio

(10)

式中:sol_sta(c,b)為校正后第b個HRU中第c層土壤的含水量,sol_st(c,b)為未校正前第b個HRU中第c層土壤的含水量,ratio為土壤水校正比率。

為避免HRU尺度上土壤含水量的過度更新,將土壤水校正比率設置在0.8～1.2之間,將土壤含水量設置在1×10-6mm至飽和土壤含水量之間。

3.3 誤差概化

降水是SWAT輸入數據中最重要的數據之一。本文假設降水輸入誤差服從對數正態分布[23],且降雨輸入誤差在時間和空間上獨立,即:

P′=Pεp(εp～lnN(μ,σ2))

(11)

式中:P′為降雨經誤差擾動后的降雨集合;P為站點實測降水量;εp為降雨誤差擾動項,服從對數正態分布,均值為μ,方差為σ2。參考以往研究設定μ=1,σ2=25%[24]。

采用比例因子法確定徑流觀測誤差協方差矩陣[25]:

Rk=(αYk)2

(12)

式中:Rk為k時刻的觀測誤差協方差矩陣;Yk為k時刻的觀測向量;α為比例因子,參考以往研究,設定α=0.1[25]。

4 結果與分析

4.1 滯后時間窗口確定

基于EnKS方法的同化模型性能受滯后時間窗口大小影響。隨著時間窗口的增大,需要考慮更多與當前狀態變量無關的誤差信息,若低估窗口長度,則損失部分滯后徑流信息。同時最佳時間窗口長度在枯水期和豐水期會發生變化。以七里街站為例,選取2018年1月1日—2018年4月1日為枯水期,2019年5月1日—2019年8月1日為豐水期,計算不同時期模型徑流量隨滯后時間窗口長度變化的模擬精度,結果見圖4。

圖4 不同時期和時間窗口長度下七里街站同化結果

由圖4可以看出,七里街站點枯水期最佳時間窗口長度為5 d,豐水期最佳時間窗口長度為1 d。豐水期最佳時間窗口長度明顯短于枯水期,這是因為豐水期河道徑流量大,SWAT模型中河道徑流流速隨水位升高而快速增加,徑流的滯后性降低,最佳時間窗口長度縮短。因此針對枯水期和豐水期的徑流數據同化,應分別選取較長和較短的時間窗口。

2018年1月1日—2020年12月31日包含枯水期和豐水期,為確定該時期最佳時間窗口,本文對比了該時期不同站點在不同時間窗口長度下的逐日徑流模擬精度,如圖5所示。

圖5 不同時間窗口長度下各水文站點同化結果

從圖5中發現模型總體精度隨時間窗口的增加呈現先增加后下降的趨勢。當同化時間窗口較小時,隨著時間窗口長度的增加,模型能夠獲取更多滯后徑流信息,從而更好地改善結果。然而當同化時間窗口過長時,后續徑流模擬誤差并不與當前水文單元土壤水誤差直接相關,主要由其他誤差導致,而模型仍將后續誤差直接歸因于當前狀態變量誤差,使得同化模型模擬精度下降。綜合豐、平、枯水文周期,最終確定七里街站點最佳時間窗口長度為2 d,沙縣站點最佳時間窗口長度為5 d,竹岐站點最佳時間窗口長度為4 d。

4.2 徑流量模擬

基于率定后的SWAT模型,利用EnKS方法和EnKF方法分別模擬閩江流域2018—2020年逐日徑流量,并與原模型模擬結果進行對比,如圖6所示,圖6(d)、6(e)、6(f)中的縱坐標為以10為底數的對數坐標。由圖6可見,原模型在枯水期徑流量模擬中存在明顯的低估現象,在豐水期則存在不同程度的低估和高估現象。而EnKS方法和EnKF方法模擬結果更為接近觀測值,表明同化站點實測徑流量能夠有效校正模型狀態誤差。

圖6 2018-2020年各水文站不同同化方法及原模型徑流量模擬結果

相較于EnKS方法,EnKF方法存在明顯的時間滯后現象。以2018年初為例,如圖7所示,EnKF方法的徑流峰值出現在觀測峰值之后,形成虛假峰值。這是由于SWAT模型中土壤水和徑流之間的映射存在滯后,在枯水期土壤水校正對當天徑流的影響極小,導致EnKF方法低估了土壤水變化的影響,對土壤水的校正不足。此外,當土壤水校正對當天徑流產生影響時,EnKF方法僅能獲取當天徑流變化信息,忽視了其余滯后徑流信息,導致EnKF方法對土壤水過度校正,出現虛假峰值。EnKS方法綜合考慮了滯后時間窗口內的模型模擬徑流和觀測徑流,能在觀測時間前對土壤含水量進行校正,模擬結果更為準確。

圖7 2018年1—4月各水文站不同同化方法及原模型徑流模擬結果

表3為不同同化方法及原模型對2018—2020年各水文站點徑流模擬精度的評價。由表3中的評價參數可知,相較原模型模擬結果,EnKS方法和EnKF方法均顯著改善了徑流模擬結果。其中EnKF方法模擬的七里街、沙縣、竹岐3個站點徑流NSE分別提升了0.06、0.03、0.01,RMSE分別降低了10.29%、5.43%、8.15%;EnKS方法模擬的七里街、沙縣、竹岐3個站點徑流NSE分別提升了0.09、0.15、0.04,RMSE分別降低了16.95%、30.78%、12.05%。3個站點NSE均達到0.80以上,其中沙縣站點效果最好,NSE達到0.86,表明EnKS方法在SWAT模型中同化站點徑流數據能夠顯著改善徑流估計精度。EnKF方法雖然同樣改善了徑流估計值,但改進程度不及EnKS方法。EnKS方法模擬結果的NSE在七里街、沙縣、竹岐3個站點上比EnKF方法分別提升了0.03、0.12、0.03,RMSE分別降低了7.43%、26.81%、4.25%。

表3 2018—2020年各水文站不同方法徑流模擬精度評價

4.3 校正土壤水對不同水文變量的影響

土壤水的校正影響地表徑流、壤中流等其他水文組分,分析土壤水校正對不同水文組分的影響,能夠為水文數據同化的狀態變量選擇提供一定參考。選擇典型HRU分析土壤水校正對其他水文變量的影響,其下墊面性質如表4所示。

表4 所選典型HRU下墊面性質

圖8為原SWAT模型和EnKS校正土壤水后模型變量基流(淺層地下水補給)、深層地下水補給、壤中流、地表徑流在2018年1月1日—2020年12月31日的變化情況。

由圖8可見,閩江流域SWAT模型中水文響應單元產流主要由基流、壤中流和地表徑流組成,深層地下水補給在產流中占比極少,且受土壤水變化影響不大,因此深層地下水補給并不是水文數據同化的理想狀態變量。而在坡度平緩地區,土壤水以下滲為主,較少產生側向徑流,壤中流在產流中的占比較少,因此對于坡度平緩地區,壤中流并不是理想的同化狀態變量。

以2018年1月1日—2018年4月1日為枯水期,2019年5月1日—2019年8月1日為豐水期,分別計算同化前后HRU不同徑流分量值在枯水期和豐水期之和,結果如表5所示。

表5 所選典型HRU在枯水期和豐水期同化前后不同徑流分量 mm

表5中的計算結果表明,SWAT模型中土壤水的校正主要改變了模型中的基流、壤中流和地表徑流。在枯水期,土壤水的校正主要改變了基流和壤中流;在豐水期,土壤水的校正主要改變了基流、壤中流和地表徑流。

對比HRU-1167、HRU-1168、HRU-1169 3種不同坡度典型下墊面,發現壤中流隨坡度的增加而增加,SWAT模型中壤中流的計算與坡度成正比,因此當流域地勢較為平緩時,土壤水校正對壤中流的改進不明顯。分別對比高坡度下墊面HRU-1166、HRU-1169和低坡度下墊面HRU-1205、HRU-1279,發現土壤性質差異同樣影響土壤水校正對壤中流的改進程度。由于SWAT模型中壤中流計算與飽和水力傳導系數SOL_K成正比,而不同類型土壤的飽和水力傳導系數存在差異,導致土壤水校正在高滲透率土壤區域對壤中流的改進更顯著。

5 討 論

EnKS方法和EnKF方法均可糾正SWAT模型土壤水誤差,改善水文變量估計結果,提高SWAT模型對流量的預測準確度。實驗結果顯示基于EnKS方法的同化效果優于EnKF方法,這和Li等[9]在GR4H水文模型以及Meng等[19]在新安江模型中所進行的實驗結果相一致,表明考慮水文模型的時間滯后性可以有效提升模型同化精度。

時間窗口的長度對EnKS方法的性能影響較大。為減少計算量,本文以比率的方式對土壤水進行校正,這需要綜合不同HRU中土壤水至河道出口的時間,以確定最佳時間窗口長度。由于SWAT模型中徑流流速與徑流量相關,因此最佳時間窗口長度隨徑流量的變化而變化。此外,由于不同流域下墊面性質的差異,最佳時間窗口在不同流域也會發生變化,這為快速準確地確定最佳時間窗口長度帶來了一定困難。本文綜合了2018年1月1日—2020年12月31日徑流數據同化總體精度,確定了閩江流域不同區域的最佳同化時間窗口長度。后續可以考慮采取自適應方法,自動確定不同時期不同區域最佳時間窗口長度,以更好地校正模型土壤含水量,改善徑流估計結果。

由于水文模型機制和下墊面性質的空間異質性,土壤水校正對不同徑流分量的改進程度在空間上存在差異。雷芳妮[26]的研究同樣表明,土壤水校正在高滲透率土壤類型上對于壤中流等水文變量能取得更顯著的改進。除土壤類型外,坡度差異也會導致土壤水校正對壤中流的改進程度不同。此外,數據同化算法對不同徑流分量的改進程度存在時間異質性,在枯水期土壤水校正主要改變了基流和壤中流,而在豐水期主要改變了基流、壤中流和地表徑流。因此,選擇地表徑流和壤中流等相關狀態變量作為同化對象時,應考慮水文周期和區域下墊面性質,以避免與模型機制相沖突。

6 結 論

(1)EnKS最優時間窗口長度在不同水文周期和不同流域存在差異。綜合豐、平、枯水文周期,確定七里街站上游最佳時間窗口長度為2 d,沙縣站上游最佳時間窗口長度為5 d,竹岐站上游最佳時間窗口長度為4 d。

(2)對比不同同化方法,EnKS方法可以有效解決水文模型的時間滯后性,減少虛假峰值現象,提高模型同化精度。該方法對徑流量模擬結果的NSE在七里街、沙縣、竹岐3個站點較EnKF方法分別提升了0.03、0.12、0.03,RMSE分別減小了7.43%、26.81%、4.25%。

(3)數據同化方法對不同徑流分量的改進程度存在空間異質性和時間異質性。在高滲透率土壤、陡坡區域,EnKS方法能使壤中流獲得更顯著的改進。在豐水期EnKS方法對地表徑流的改進較枯水期更為明顯。