通用線性模型在氣象水文集合預報后處理中的應用

姜迪，智海，趙琳娜，段青云，梁莉，劉瑩

(1．南京信息工程大學，江蘇南京210044;2．中國氣象科學研究院災害天氣國家重點實驗室，北京100081;3．北京師范大學，北京100875;4．中國氣象局 公共氣象服務中心，北京100081;5．四川省氣象臺，四川 成都610072)

0 引言

洪水災害是世界上危害最大、造成損失最多的自然災害之一。但是由于其過程復雜，包含的不確定因素較多，目前的確定性預報對降水(氣象)和流量(水文)的預報都難以滿足實際需要，因此概率預報不失為一條提高水文預報精度的途徑。采用概率預報的形式對水文不確定性進行定量描述，不僅可以提高預報的精度和可靠性，而且可以為決策者做出客觀的最優(yōu)決策提供更豐富的信息，同時也能滿足用戶對風險信息的需求。隨著集合概率預報得到越來越廣泛的使用，它有可能會取代單一預報(杜鈞，2002;智協(xié)飛和陳雯，2010)。

本文的后處理是指在氣象—水文集合預報系統(tǒng)中消除水文模型模擬過程產(chǎn)生的誤差，降低水文模擬中不確定性影響的過程。水文后處理方面的研究以貝葉斯方法取得的成果較為突出。Krzysztofowicz(1985)首先提出了基于正態(tài)線性假設的貝葉斯預報處理器(Bayesian Processor of Forecasts，簡稱BPF)，通過貝葉斯公式求得待預報流量的后驗分布，并基于后驗分布提出最優(yōu)決策。隨后，Krzysztofowicz(1999)又提出了基于貝葉斯理論框架的水文預報系統(tǒng)(Bayesian Forecasting System，簡稱BFS)，將降水預報不確定性與水文模型不確定性分別量化，然后通過全概率公式將二者耦合起來，得到水文預報不確定性的解析解，BFS可以與任何模型或預報方案相結合使用。之后Krzysztofowicz and Kelly(2000)在BPF的基礎上提出了水文不確定性處理器(Hydrologic Uncertainty Processor，簡稱 HUP)，HUP的核心是亞高斯模型，對實測值和預報流量進行正態(tài)分位數(shù)轉化(使得轉化變量的線性正態(tài)關系更加顯著)，由貝葉斯公式得到轉化變量的后驗密度函數(shù)，然后還原得到待預報流量的后驗密度函數(shù)。Krzysztofowicz and Maranzano(2004)還提出了一種基于HUP方法產(chǎn)生的概率水位預報來消除誤差?；谶@種方法，Seo et al．(2006)為美國國家氣象局的河流預報系統(tǒng)提供了一種水文后處理器。國內對水文預報的不確定性也做了很多探索，王善序(2001)、錢明開等(2004)、張洪剛和郭生練(2004)、張洪剛等(2004)、梁忠民等(2010)利用 BFS、HUP等方法進行了水文不確定性和概率水文預報的研究，結果表明，概率預報與確定性預報同樣有價值，特別是當預報不確定性比較大時，概率預報比確定性預報在決策方面具有更高的經(jīng)濟價值。貝葉斯方法在氣象集合預報中也有諸多應用，陳法敬等(2011)及段明鏗和王盤興(2006)利用貝葉斯方法提高了集合預報系統(tǒng)的性能，提高了預報精度。另外，還有一些其他研究方法，如:Georgakakos and Smith(1990)利用卡爾曼濾波技術與水文模型相結合，得到模擬流量的均值與方差，進行短期洪水預報;Georgakakos et al．(1998)通過蒙特卡洛法隨機生成若干組不同代表類型的氣象資料，以不同的初始條件輸入水文氣象模型中得到所有可能的徑流過程，延長了預見期，提高了模擬的可見性。

但是上述方法都有各自的缺陷，例如貝葉斯方法需要知道不確定性的概率分布和協(xié)方差結構，卡爾曼濾波方法無法輸出不確定性的解析解，蒙特卡洛方法在隨機抽取樣本方面收斂速度慢，精確度不高，對于三維以下的問題不如基于線性最小二乘回歸的自適應方法、分對數(shù)回歸模型的自適應算法和卡爾曼濾波方法好(馮圓和龔曉燕，2010)。而通用線性模型(General Linear Model，簡稱GLM)方法可以避開這些問題，得到比較可靠的概率水文預報。

1 通用線性模型方法

從概率上看，水文后處理過程實際上是要得到一個條件概率密度函數(shù)，這個條件概率密度函數(shù)是f(yobs|yfcst)，即給定了氣象預報yfcst的觀測值yobs條件概率密度函數(shù)。如果忽略其他原因產(chǎn)生的誤差，

只關注于模型模擬水文事件產(chǎn)生的誤差，通用線性模型后處理方法的原理可以用方程(1)表示:

圖1是一個GLM的數(shù)據(jù)窗口示意圖。根據(jù)Zhao et al．(2011)對通用線性模型方法的研究，數(shù)據(jù)窗口主要分為分析期和預報期。分析期的作用是提取資料樣本以確定GLM模型參數(shù)，然后使用確定好的GLM模型對預報期的資料進行校正，得到校正過的預報值。Na和Nf分別代表分析期和預報期的長度。在分析期和預報期的前后各有一個buffer期，作用是增加參與確定GLM模型參數(shù)的樣本數(shù)量，更好地確定 GLM模型參數(shù)。Nbuffer表示buffer期的長度。

如果隨機變量yobs(觀測值)、ysim(模擬值)和yfcst(預報值)是服從正態(tài)分布的，就可以利用通用線性模型方法去解決方程(1)的問題(R and Schakke，1973)。

圖1 通用線性模型的數(shù)據(jù)窗口示意圖Fig．1 Sketch map of data window of GLM

由于上述方法要求Z1和Z2中的變量服從正態(tài)分布，顯然觀測和模擬流量并不滿足這一條件，所以本文使用由Krzysztofowicz(1985)提出的正態(tài)分位數(shù)轉換方法(Normal Quantile Transformation，簡稱NQT)來解決這一問題。這種方法是將原始空間中的流量數(shù)據(jù)轉換為在轉換空間中的正態(tài)分布變量，使轉化變量的線性正態(tài)關系更加顯著。隨后進行一個逆NQT過程，將轉換變量的估計值轉換回原始空間才能使用并符合正態(tài)分布的數(shù)據(jù)。

2 數(shù)據(jù)與試驗介紹

本文使用的數(shù)據(jù)來自國際模型參數(shù)估計試驗(ModelParameterEstimationExperiment，簡 稱MOPEX)(Duan et al．，2005)。MOPEX 是一個關于水文模型和大氣模式的陸面參數(shù)化方案的國際合作項目，得到了世界氣象組織(WMO)、國際水文科學協(xié)會(簡稱IAHS)等多個氣象、水文部門和科學組織的支持。它通過利用氣象資料和水文資料對大量水文模型在相當數(shù)量的河流流域里進行模型參數(shù)的率定試驗，以探索模型最優(yōu)參數(shù)的估計方法。入選MOPEX的流域無大型水利工程，受人工影響較少，可以反映天然狀況下的水文過程。美國本土與中國緯度相近，有多種相似的氣候類型，尤其是密西西比河流域與中國長江流域有很多共同點，可以為中國的水文研究提供參考。

MOPEX中水文模型輸出的模擬數(shù)據(jù)分為兩種，一種為經(jīng)過模型參數(shù)率定的水文模型產(chǎn)生的模擬數(shù)據(jù)(簡稱CAL數(shù)據(jù))，其模擬誤差較小;另一種為未經(jīng)過參數(shù)率定的水文模型產(chǎn)生的模擬數(shù)據(jù)(簡稱APR數(shù)據(jù))，其模擬誤差較大。

本文選取MOPEX數(shù)據(jù)庫中French Broad River流域作為試驗流域。此流域面積為2 448 km2，此區(qū)域73%的覆蓋面積為森林，年降水量為1 676 mm。同時選取薩克拉門托水文模型(Sacramento，簡稱SAC模型)作為試驗模型，數(shù)據(jù)的時間長度為36 a，即1962年1月1日至1997年12月31日。SAC模型最早是一個設計用于洪水預報的概念性的集總式模型。它包括16個參數(shù)，利用一系列有一定物理概念的數(shù)學公式來描述水分的運動過程，具有較強的物理概念和廣泛的適用性(劉金平和樂嘉祥，1996)。目前SAC模型在美國被廣泛應用于水文預報，同時也是中國國內引入較早、研究應用比較多的水文模型。

確定了模型和流域后，本文將根據(jù)不同的試驗設計選取分析期資料確定GLM模型參數(shù)，然后用確定好參數(shù)的GLM模型去校正預報期的流量模擬值并產(chǎn)生關于預報期的集合預報，依據(jù)對不同誤差水平的兩種資料(APR資料和CAL資料)的后處理效果來檢驗GLM校正誤差的能力:1)使用GLM分別對兩種模擬數(shù)據(jù)進行校正，分析期和預報期均與資料時間長度相同，通過若干統(tǒng)計指標將校正結果與實際觀測對比，檢驗誤差校正效果;2)GLM使用不同的校驗策略，即使用不同年代長度的分析期資料來確定GLM模型參數(shù)，然后對預報期的模擬資料進行后處理，測試在使用不同數(shù)量樣本確定的模型參數(shù)下GLM校正誤差的能力。

3 結果分析

3.1 兩種數(shù)據(jù)集GLM在不同誤差水平下的誤差校正能力

試驗中，設定數(shù)據(jù)窗的參數(shù)為:Na取30 d，Nf取30 d，Nbuffer取30 d;預報起始日期為2月24日。試驗選用36 a資料確定GLM模型參數(shù)，并用36 a資料進行驗證。

首先，根據(jù)第2節(jié)所述，需對數(shù)據(jù)進行NQT轉換，使其符合正態(tài)分布。圖2為APR數(shù)據(jù)和CAL數(shù)據(jù)在NQT轉換前后與實際觀測的相關系數(shù)比較?？梢?，在進行NQT轉換前，CAL數(shù)據(jù)與實際觀測的相關系數(shù)基本達0.8，而APR數(shù)據(jù)與實際觀測的相關系數(shù)為0.6～0.8，CAL數(shù)據(jù)的相關性明顯好于APR數(shù)據(jù)。經(jīng)NQT轉換后，兩數(shù)據(jù)與實際觀測的相關系數(shù)均有提高，CAL數(shù)據(jù)的相關系數(shù)為0.8～0.9(達0.9的數(shù)據(jù)明顯較多)，仍然好于APR數(shù)據(jù)(僅個別數(shù)據(jù)的相關系數(shù)達0.9)。

圖3是使用CAL資料和APR資料產(chǎn)生的集合流量預報的連續(xù)等級概率評分(Continuous Rank Probability Score，簡稱CRPS)的比較。CRPS的值可以定量地評價集合預報的性能;值越大表示集合預報的性能越差，值越小表示集合預報越好，零分代表最好預報。由圖3可見，兩種資料產(chǎn)生的集合流量預報的CRPS值在整個預報期內都處于較低水平，使用APR資料的集合預報的CRPS值全部低于0.45，而使用CAL資料的集合預報的CRPS值都低于0.4，說明GLM后處理器產(chǎn)生的集合流量預報是較可靠的。使用兩種資料的集合預報在預報期的初始階段均較好，得到了很低的CRPS值，隨著預報時效延長，CRPS值逐漸升高，說明預報效果不斷變差，這也是符合預報的一般規(guī)律的。整體而言，可以認為整個預報期內的集合預報效果是可靠的。

圖2 NQT轉換前后原始預報與實際觀測的相關系數(shù) a．APR資料轉換前;b．APR資料轉換后;c．CAL資料轉換前;d．CAL資料轉換后Fig．2 Correlation coefficients between raw simulations and observations(before or after NQT transformation) a．APR data before NQT transformation;b．APR data after NQT transformation;c．CAL data before NQT transformation;d．CAL after NQT transformation

圖3 使用APR和CAL資料的集合預報的CRPS比較Fig．3 Comparison of CRPS scores from ensemble forecasts using APR and CAL data

圖4是GLM使用兩種資料的校正結果與實際觀測在統(tǒng)計指標上的對比。在下面的對比中選用集合預報結果的集合平均作為單值的校正預報，這是因為集合平均位于集合中位與預報分布的中心，代表了真實情況下最有可能出現(xiàn)的狀況。

圖4 使用不同資料水文模擬結果的評估 a．原始模擬、校正值與實際觀測的平均值對比;b．原始模擬、校正值與實際觀測的標準差對比;c．原始模擬和校正值的均方根誤差對比;d．APR資料的原始模擬、校正值與實際觀測的累積分布函數(shù)比較;e．CAL資料的原始模擬、校正值與實際觀測的累積分布函數(shù)比較Fig．4 Evaluations of results using different simulation data a．comparison of mean of raw simulation，adjusted and observation;b．comparison of STD of raw simulation，adjusted and observation;c．comparison of RMSE of raw simulation and adjusted;d．comparison of CDF of raw simulation，adjusted and observation using APR data;e．comparison of CDF of raw simulation，adjusted and observation using CAL data

圖4a是使用CAL資料和APR資料的原始模擬值、校正預報值與實際觀測三者平均值(簡稱Mean)的比較。校正預報的流量是否與實際觀測相吻合是對校正效果好壞最直觀的判斷。圖4a表明，不論使用APR資料還是使用CAL資料，預報期內校正預報的平均都比原始模擬的平均更加接近于實際觀測平均。這表明經(jīng)過GLM后處理器的處理，有誤差的水文模擬得到了很好的校正，結果更接近實際觀測。對APR資料和CAL資料的校正結果進行對比可知，誤差較小的CAL資料的校正結果相對于誤差較大的APR資料的校正結果更接近實際觀測，說明GLM不僅可以處理誤差較大的水文模擬，對于誤差較小的水文模擬也有一定的處理能力。

圖4b是使用CAL資料和APR資料的原始模擬值、校正預報值與實際觀測三者的標準差(Standard Deviation，簡稱STD)的比較。STD是對不確定性的一種量化。STD越高，表示試驗數(shù)據(jù)越離散，即觀測值與預測值相差越大;反之，STD越低，則代表觀測值與預測值相差越小。圖4b表明，不論是使用APR資料還是使用CAL資料的結果，經(jīng)過GLM后處理的校正預報的標準差都比原始模擬更接近于實際觀測，尤其是對APR資料的改善極其顯著。其中盡管APR資料的校正結果的STD相比于CAL資料的校正結果的STD得到了較大程度的校正，更接近于實際觀測，但由于被校正得過低，反而與真實情況有差距。對比CAL資料的標準差的校正情況，可以認為APR資料經(jīng)過GLM校正后，還包含了較大的不確定性，這些不確定性可能與水文模型自身的參數(shù)率定有關。

圖4c是使用CAL資料和APR資料的校正預報值與原始模擬值的均方根誤差(Root-Mean-Square Error，簡稱RMSE)的比較。RMSE是由實際觀測值與預報值之間的差值定義的，RMSE的值越大，表示預報誤差越大，反之則表明預報誤差越小?？梢钥吹?，相對于原始預報的RMSE，經(jīng)過后處理的校正預報的RMSE明顯降低。使用APR資料的校正預報的誤差被降低到使用CAL資料的原始預報的誤差幾乎同樣水平，而CAL資料的校正預報的誤差則得到了更進一步降低，這也驗證了之前分析的結論。

圖4d和圖4e分別是APR資料和CAL資料的原始模擬、校正預報與實際觀測的累積分布函數(shù)的比較。累積分布函數(shù)(Cumulative Distribution Function，簡稱CDF)描述的實數(shù)隨機變量x的概率分布，也是概率密度函數(shù)(Probability Density Function，簡稱PDF)的積分。概率密度函數(shù)表示的是隨機變量的瞬時取值落在某指定范圍內的概率，CDF代表的是隨機變量x小于或者等于某個數(shù)值的概率，也就是連續(xù)函數(shù)F(x)的分布。圖4d和圖4e表明，APR資料和CAL資料經(jīng)過GLM處理的校正值在分布上都更接近實際觀測;特別是CAL資料，經(jīng)過水文模型自身的參數(shù)率定之后，其誤差已得到較大程度的降低，數(shù)據(jù)分布與實際觀測數(shù)據(jù)相當接近，再經(jīng)過GLM后處理后，數(shù)據(jù)分布狀況得到進一步改善。

綜合所述，這些試驗初步證明GLM不僅可以有效降低水文模擬中的誤差，還可以產(chǎn)生可靠的集合流量預報。它不僅對誤差較大的水文模擬能取得較好的校正效果，而且對誤差較小的水文模擬也有一定的處理能力，可以進一步降低誤差。

3.2 不同的校驗策略下GLM的性能比較

在下面的試驗中，數(shù)據(jù)窗口的參數(shù)不變，只是選取了另外兩種不同的校驗GLM模型參數(shù)的策略(即選擇不同年代長度的分析期和驗證期資料)，來進行與第3.1節(jié)一樣的試驗，以測試GLM在不同校驗策略下的校正效果。

由表1可見，校驗期策略1(即第3.1節(jié)使用的試驗設計)表示用36 a長度的資料確定GLM模型參數(shù)，然后用同樣的36 a資料進行驗證;校驗期策略2表示用前18 a資料確定GLM模型參數(shù)，用后18 a資料進行驗證;校驗期策略3表示用單數(shù)年(共18 a)資料確定GLM模型參數(shù)，用偶數(shù)年(共18 a)資料進行驗證。

圖5a、圖6a分別是在使用校驗策略2、3情況下用不同資料得到的模擬值、校正值與實際觀測平均值的比較。可見，不管使用哪種校驗策略，經(jīng)過GLM處理的校正值都比原始模擬值要更加接近實際觀測，誤差得到有效降低。而且使用CAL資料的校正值在兩個統(tǒng)計指標上也都比使用APR資料的校正值更加接近實際觀測，其他指標如標準差、累積分布函數(shù)等的結果也大致相同，不再贅述。該結果既支持了上一個試驗所得到的結論，又說明GLM后處理器在不同校驗策略下(即在不同的參數(shù)配置下)，對水文模擬中不確定性的校正作用也是有效的。

圖5 校驗策略2的結果 a．原始模擬、校正值與實際觀測的平均值對比;b．原始模擬和校正值的均方根誤差對比Fig．5 Results of the calibration strategy 2 a．comparison of mean of raw simulation，adjusted and observation;b．comparison of RMSE of raw simulation and adjusted

圖6 校驗策略3的結果 a．原始模擬、校正值與實際觀測的平均值對比;b．原始模擬和校正值的均方根誤差對比Fig．6 Results of the calibration strategy 3 a．comparison of mean of raw simulation，adjusted and observation;b．comparison of RMSE of raw simulation and adjusted

但是，也應注意到，在預報的某些時段內，使用CAL資料的校正預報誤差并不比使用APR資料的校正預報誤差更低(圖5b、圖6b)，這表明并不是所有的有誤差的水文模擬通過GLM后處理都能得到改善。其原因在于:如果對降水資料的前處理以及水文模型自身的參數(shù)率定等前期工作可以得到比較理想的結果，使得水文模擬中的誤差已經(jīng)消除得比較多，則后處理過程無法再進一步減小誤差。

由于3種校驗策略的本質區(qū)別是參與確定GLM模型參數(shù)的樣本數(shù)量不同，可以看到使用策略1的樣本數(shù)量最多，使得使用策略1得到的校正效果最好，所以參與確定GLM模型參數(shù)的樣本值越多，提取的歷史信息越多，校正的結果就越接近實際觀測。

4 結論

本文利用MOPEX提供的兩種不同誤差水平的資料，對比檢驗GLM后處理器對水文集合預報誤差的校正能力。結果如下:

1)GLM后處理器可以產(chǎn)生可靠的水文集合預報，集合預報的連續(xù)分級概率評分值都在比較低的水平，證明預報結果是可靠的。

2)GLM后處理器對有偏差的水文預報的校正效果是顯著的。在平均值、標準差、均方根誤差等統(tǒng)計指標方面的比較顯示，無論是誤差較大或者是誤差較小的水文模擬，校正預報都要比原始模擬更加接近于實際觀測，誤差顯著降低。

3)即使是對經(jīng)過參數(shù)率定的水文模型產(chǎn)生的誤差較小的水文模擬，GLM后處理器仍然有能力消除一部分誤差。

4)使用不同的校驗策略，選取不同的分析期資料對不同預報期資料進行校正，GLM后處理器都能很好地反映出實際觀測的變化趨勢，降低原始模擬中的誤差。這表明GLM的處理效果是穩(wěn)定有效的，不會因為較少的資料而造成誤差校正效果的顯著下降。

5)通過試驗結果對比，發(fā)現(xiàn)當原始模擬已經(jīng)做得足夠好，模擬誤差已經(jīng)比較小的時候，后處理不一定能得到更好的結果。

本文只選取了MOPEX數(shù)據(jù)庫中一個模型在一個流域上的資料進行試驗，初步驗證了GLM的性能。下一步要擴大流域和模型的數(shù)量，特別是加強國內流域和水文模型的試驗，以進一步證明GLM的校正能力是普遍有效的。

陳法敬，矯梅燕，陳靜．2011．一種溫度集合預報產(chǎn)品釋用方法的初步研究［J］．氣象，37(1):14-20．

杜鈞．2002．集合預報的現(xiàn)狀和前景［J］．應用氣象學報，13(1):16-28．

段明鏗，王盤興．2006．一種新的集合預報權重平均方法［J］．應用氣象學報，17(4):488-493．

馮圓，龔曉燕．2010．基于蒙特卡羅方法的氣象問題應用研究［C］//第27屆中國氣象學會年會大氣物理學與大氣環(huán)境分會場論文集．北京:中國氣象學會．

梁忠民，戴榮，王軍，等．2010．基于貝葉斯模型平均理論的水文模型合成預報研究［J］．水力發(fā)電學報，29(2):114-118．

劉金平，樂嘉祥．1996．薩克拉門托模型參數(shù)初值分析方法研究［J］．水科學進展，3(7):252-259．

錢名開，徐時進，王善序，等．2004．淮河息縣站概率預報模型研究［J］．水文，24(2):23-25．

王善序．2001．貝葉斯概率水文預報簡介［J］．水文，21(5):33-34．

張洪剛，郭生練．2004．貝葉斯概率洪水預報系統(tǒng)［J］．科學技術與工程，4(2):74-75．

張洪剛，郭生練，劉攀，等．2004．基于貝葉斯分析的概率洪水預報模型研究［J］．水電能源科學，22(1):22-25．

智協(xié)飛，陳雯．2010．THORPEX國際科學研究新進展［J］．大氣科學學報，33(4):504-511．

Duan Q，Schaake J，Andreassian V，et al．2005．Model parameter estimation experiment(MOPEX):Overview and summary of the second and third workshop results［J］．J Hydrol，320(1/2):3-17．

Georgakaos K P，Smith G F．1990．On improved hydrologic forecasting—Results from a WMO real-time forecasting experiment［J］．J Hydrol，114(1):17-45．

Georgakakos A P，Yao H，Mullusky M G，et al．1998．Impacts of climate variability on the operational forecast and management of the upper Des Moines River basin［J］．Water Resour Res，34(4):799-821．

Krzysztofowicz R．1985．Bayesian model of forecasted time series［J］．Water Resour Res，21(5):805-814

Krzysztofowicz R．1999．Bayesian theory of probabilistic forecasting via deterministic hydrologic model［J］．Water Resour Res，35(9):2739-2750．

Krzysztofowicz R，Kelly K S．2000．Hydrologic uncertainty processor for probabilistic river stage forecasting［J］．Water Resour Res，36(11):3265-3277．

Krzysztofowicz R，Maranzano C J．2004．Hydrologic uncertainty processor for probabilistic stage transition forecasting［J］．J Hydrol，293(1):57-73．

Seo D J，Herr H D，Schaake J C．2006．A statistical post-processor for accounting of hydrologic uncertainty in short-range ensemble streamflow prediction［J］．Hydrology and Earth System Sciences Discussions，3(4):1987-2035．

R D V，Schakke J C Jr．1973．Disaggregation processes in stochastic hydrology［J］．Water Resour Res，9(3):580-585．

Zhao L，Duan Q，Schaake J，et al．2011．A hydrologic post-processor for ensemble streamflow predictions［J］．Adv Geosci，29:51-59．