廣義相加模型在烏江夏季徑流預(yù)報(bào)中的應(yīng)用

榮艷淑，胡玉恒，馮瑞瑞，殷雨婷，李崇浩

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098；2.河海大學(xué)水利學(xué)科專業(yè)實(shí)驗(yàn)教學(xué)中心，江蘇 南京 210098；3. 南京恩瑞特實(shí)業(yè)有限公司，江蘇 南京 211106；4.中國南方電網(wǎng)電力調(diào)度控制中心，廣東 廣州 510623)

中長期徑流預(yù)測(cè)一直是水文預(yù)報(bào)的一大難題。傳統(tǒng)的統(tǒng)計(jì)水文預(yù)報(bào)模型，如時(shí)間序列分析和各種回歸模型等，由于計(jì)算簡(jiǎn)便，模型參數(shù)有一定的物理意義而被廣泛地應(yīng)用于水文預(yù)報(bào)工作[1]。如，基于多種氣象因子建立的逐步回歸模型、多元線性回歸模型和時(shí)間序列分解模型對(duì)月、旬徑流都有較好的預(yù)報(bào)效果[2-4]。但是傳統(tǒng)線性統(tǒng)計(jì)模型沒有考慮地球氣候系統(tǒng)中各要素與徑流之間復(fù)雜的非線性關(guān)系，預(yù)測(cè)精度很難有進(jìn)一步提高[5]。隨著計(jì)算機(jī)技術(shù)的發(fā)展，人工神經(jīng)網(wǎng)絡(luò)、改進(jìn)人工神經(jīng)網(wǎng)絡(luò)和小波分解等，不僅能夠處理變量之間的非線性映射關(guān)系[6]，而且預(yù)測(cè)結(jié)果明顯優(yōu)于線性回歸模型，預(yù)報(bào)精度也有明顯提高[7-9]。但是這些方法存在模型物理意義不明確，模型通用性較差等缺陷[10]。近年來涌現(xiàn)了一些考慮了非線性關(guān)系，或者使用多種模型融合的技術(shù)方法，明顯提高了預(yù)測(cè)水平[11-15]。廣義相加模型(generalized additive model,GAM)就是綜合了線性關(guān)系和非線性關(guān)系的一種復(fù)合模型，它承認(rèn)預(yù)測(cè)變量與某些物理因子存在線性關(guān)系，與另外的物理因子可能存在非線性關(guān)系，并基于線性和非線性因子建立綜合模型。這一模型在基于海溫因子的降水預(yù)測(cè)、基于氣溫和降水的大氣污染影響分析和基于NAO的阿爾卑斯山冬季極端氣溫的預(yù)測(cè)中都得到較好的應(yīng)用[16-18]。

烏江流域位于長江上游，流域內(nèi)水電站網(wǎng)密布，為當(dāng)?shù)氐霓r(nóng)業(yè)灌溉、生活用水以及防洪發(fā)電帶來了顯著效益。但是，流域內(nèi)徑流預(yù)測(cè)準(zhǔn)確性一直困擾著水電站網(wǎng)管理部門，特別是汛期季風(fēng)盛行期間，降水變化顯著，徑流變化更難以預(yù)測(cè)。在這一背景下，基于廣義相加模型，筆者對(duì)烏江流域的洪家渡夏季徑流進(jìn)行預(yù)測(cè)研究，并對(duì)GAM的建立、模擬和預(yù)測(cè)進(jìn)行評(píng)估和探討，以期能為烏江流域的徑流預(yù)測(cè)提供一個(gè)新途徑。

1 研究數(shù)據(jù)與方法

1.1 研究數(shù)據(jù)

以烏江流域的洪家渡夏季徑流為研究對(duì)象，以多種海溫指數(shù)作為預(yù)測(cè)因子，包括太平洋、大西洋和印度洋的26個(gè)海溫指數(shù)等(表1)。徑流資料來自洪家渡水電站，海溫指數(shù)來自國家氣候中心。由于海溫與我國夏季降水之間存在較長的滯后相關(guān)性[19]，因此選取前期冬季海溫指數(shù)進(jìn)行研究。冬季指數(shù)為12月、翌年1—2月3個(gè)月的平均值，夏季徑流為6—8月3個(gè)月的平均值。資料時(shí)間范圍均為1952—2014年。

表1 海溫指數(shù)與因子代碼

1.2 GAM模型

GAM是由Hastie等[20]在廣義線性模型(GLM)[21]的基礎(chǔ)上改進(jìn)而來，該模型包括兩部分，分別描述了自變量與因變量的線性和非線性關(guān)系，因此GAM存在參數(shù)和非參數(shù)兩部分，屬于半?yún)?shù)模型。這一改進(jìn)打破了線性模型中關(guān)于自變量與因變量線性相關(guān)的假定，使得模型應(yīng)用有了進(jìn)一步拓展。GAM可表示為

(1)

式中：g(yi)——因變量yi的連接函數(shù)；l——序列長度；β0——截距；s(xki)——平滑函數(shù)，本文采用薄板樣條函數(shù)作為平滑函數(shù)；n——非線性因子的個(gè)數(shù)；xki——非線性因子；βj——線性部分的回歸系數(shù)；m——線性因子的個(gè)數(shù)；xji——線性因子；εi——隨機(jī)誤差，服從高斯正態(tài)分布。

取l= 63，為避免平滑過程中出現(xiàn)過度擬合，將薄板樣條函數(shù)節(jié)點(diǎn)的最大個(gè)數(shù)固定為5。GAM參數(shù)部分檢驗(yàn)使用t檢驗(yàn)，非參數(shù)部分使用F檢驗(yàn)。

1.3 模型模擬和預(yù)測(cè)的評(píng)估方法

1.3.1 模型模擬的評(píng)估方法

基于式(1)對(duì)洪家渡夏季徑流序列進(jìn)行模擬，模擬效果選用最小信息準(zhǔn)則(Akaike information criterion, AIC)、均方根誤差(root mean square error, RMSE)、平均絕對(duì)誤差(mean absolute error, MAE)、概率空間線性誤差(linear error in probability space, LEPS)[22]和線性相關(guān)系數(shù)(correlation,r2)5種評(píng)價(jià)指標(biāo)進(jìn)行評(píng)估。根據(jù)最小信息準(zhǔn)則，AIC值越小表明模型模擬效果越好；RMSE用來衡量模擬值與觀測(cè)值之間的相對(duì)偏差，其值越小表明偏差越小；MAE能夠反映模擬值與實(shí)測(cè)值的絕對(duì)偏離程度，其值越小模擬效果越好；LEPS使用了Potts等[23]改進(jìn)的LEPS，其值越接近于1，模擬效果越好。r2刻畫2個(gè)序列之間的線性相關(guān)，其值越接近1表明整個(gè)序列模擬值與實(shí)測(cè)值的相關(guān)性越高。

1.3.2 模型預(yù)測(cè)的評(píng)估方法

一般的模型都是取序列長度的一部分?jǐn)?shù)據(jù)建立模型，剩余數(shù)據(jù)作為預(yù)測(cè)樣本，用以評(píng)估模型的模擬能力[24]，但是這種方法存在一定的主觀性以及測(cè)試樣本長度的局限性，特別是當(dāng)預(yù)測(cè)樣本過長時(shí)，模型易出現(xiàn)衰退，導(dǎo)致預(yù)測(cè)效果不好。

交叉驗(yàn)證是一個(gè)可以克服上述問題的模型預(yù)測(cè)檢驗(yàn)方法，其主要思想是將數(shù)據(jù)樣本分成兩部分，一部分作為訓(xùn)練集進(jìn)行模型訓(xùn)練，另一部分作為測(cè)試集進(jìn)行模型預(yù)測(cè)，但是該方法是將數(shù)據(jù)循環(huán)分組，對(duì)每個(gè)原數(shù)據(jù)都進(jìn)行模擬，這樣通過多次重復(fù)試驗(yàn)?zāi)艿玫捷^好的模型估計(jì)結(jié)果[23]。為了充分利用已有的樣本數(shù)據(jù)得到更可靠的預(yù)測(cè)結(jié)果，本文使用了廣義交叉驗(yàn)證[25]確定薄板樣條函數(shù)平滑自由度的取值，使用留一法交叉驗(yàn)證[26]進(jìn)行模型預(yù)測(cè)結(jié)果的評(píng)估。

留一法交叉驗(yàn)證的主要思想是從l個(gè)樣本數(shù)據(jù)中取一個(gè)作為測(cè)試集，剩下的l-1個(gè)作為訓(xùn)練集，例如，對(duì)于某時(shí)間序列的第一個(gè)樣本，用i= 2，3，…，l個(gè)樣本建立模型，獲得模型參數(shù)，再對(duì)第一個(gè)樣本進(jìn)行模擬；對(duì)第二個(gè)樣本，則是用i= 1和i= 3，4，…，l個(gè)樣本建立模型，再對(duì)第二個(gè)樣本進(jìn)行模擬，如此往復(fù)，共進(jìn)行l(wèi)次訓(xùn)練，得到l個(gè)預(yù)測(cè)結(jié)果。最后，通過模擬序列和實(shí)測(cè)序列的相關(guān)性分析和預(yù)報(bào)相對(duì)誤差頻數(shù)統(tǒng)計(jì)，對(duì)GAM預(yù)報(bào)效果進(jìn)行評(píng)價(jià)。

表2 GAM中各預(yù)報(bào)因子的線性和非線性屬性

2 結(jié)果與討論

2.1 模型因子挑選

本文選擇前期冬季作為海溫影響的關(guān)鍵時(shí)期，將前期冬季26個(gè)海溫指數(shù)作為洪家渡夏季徑流預(yù)測(cè)的潛在預(yù)報(bào)因子。同時(shí)，為了驗(yàn)證GAM模擬的效果，使用GLM進(jìn)行對(duì)比分析。根據(jù)逐步回歸挑選因子的最小信息準(zhǔn)則[27](AIC)分別篩選出GAM和GLM模型的預(yù)報(bào)因子。GAM挑選出13個(gè)預(yù)報(bào)因子，分別為x1、x2、x3、 x4、 x6、 x7、 x10、 x14、 x17、 x20、 x22、 x24、 x26；GLM挑選出14個(gè)預(yù)報(bào)因子，分別為x4、 x5、 x7、 x8、 x9、 x10、x11、 x12、 x13、x14、 x17、 x20、x22、 x24。2個(gè)模型共有的預(yù)報(bào)因子是8個(gè)，這是線性模型和非線性模型對(duì)自變量影響因變量的物理機(jī)制存在差異的原因。

用薄板樣條函數(shù)分別對(duì)GAM中的13個(gè)因子進(jìn)行平滑處理，平滑自由度(用edf表示)的值為1時(shí)，表明自變量與因變量為線性關(guān)系，超過1時(shí)，表明自變量與因變量為非線性關(guān)系，其值越大，表明非線性關(guān)系越明顯。經(jīng)過平滑處理，有9個(gè)因子與洪家渡徑流是線性關(guān)系，有4個(gè)因子與洪家渡徑流是非線性關(guān)系(表2)。由表2可見，NINO1+2區(qū)海表溫度距平指數(shù)(x1)、黑潮海溫指數(shù)(x20)、熱帶印度洋全區(qū)一致海溫指數(shù)(x24)和副熱帶南印度洋偶極子指數(shù)(x26)均對(duì)洪家渡徑流有非線性影響，而其他指數(shù)，包括NINO3區(qū)(x2)、NINO4區(qū)(x3)、NINO3.4區(qū)(x4)、NINO C區(qū)(x6)和NINO A區(qū)(x7)海表溫度距平指數(shù)、熱帶北大西洋海溫指數(shù)(x10)、印度洋暖池強(qiáng)度指數(shù)(x14)、大西洋多年代振蕩指數(shù)(x17)和暖池型ENSO指數(shù)(x22)與洪家渡徑流的關(guān)系可用線性關(guān)系描述。

圖1 x20和x24與洪家渡夏季徑流偏差的相關(guān)關(guān)系Fig.1 Relations between summer runoff deviation at Hongjiadu and factors, x20 and for x24

以x20和x24為例，圖1給出它們與模擬徑流量偏差的相關(guān)關(guān)系，圖中實(shí)線為平滑函數(shù)值，點(diǎn)線之間的區(qū)域?yàn)樾哦?5%水平的范圍?？梢钥闯?，x20和x24與模擬值偏差的散點(diǎn)基本分布在信度95%水平之間，表明對(duì)x20和x24作為非線性因子的判定是合理的。

2.2 模型模擬效果評(píng)估

由圖2可以看出GAM和GLM對(duì)徑流量模擬的趨勢(shì)相似，且與實(shí)測(cè)序列基本相一致。從一些異常年份看，比如1954年、1957年、1967年、1998年、2001年、2008年、2012年等幾個(gè)流量異常偏多的年份，以及2005年之后幾個(gè)流量異常偏少的年份，均得到了較為準(zhǔn)確的模擬。通過GAM與GLM的對(duì)比分析，在1954年、1957年、1998年、2001年、2008年等幾個(gè)流量異常偏多的年份，GAM的模擬效果更接近實(shí)測(cè)值。在2個(gè)模型模擬值明顯不一致的年份，如1967年、1985年、1994年，GAM的模擬值與實(shí)測(cè)值更為接近，模擬效果更好。

圖2 GAM和GLM模擬值與實(shí)測(cè)值比較Fig.2 Comparison between simulated values of GAM and GLM and observed values

由表3可知， GAM的AIC值為733.04，而GLM的AIC值為752.57，表明GAM的模擬擬合效果比GLM好。GAM的RMSE與MAE值均低于GLM，說明GAM在整體序列上的模擬精度高于GLM。GAM模擬值與實(shí)測(cè)值序列的相關(guān)性(r2)與GLM相比也有很大的提升。綜合分析，GAM模擬效果比GLM模擬效果有了顯著提升。

表3 GAM和GLM評(píng)價(jià)指標(biāo)對(duì)比

由圖3(a)可以看出，模型殘差的分布大致呈正態(tài)分布的規(guī)律，而當(dāng)值大于100以及小于-50時(shí)，部分散點(diǎn)并不完全處于1∶1直線上；圖3(b)也顯示出殘差分布出現(xiàn)了略呈偏態(tài)的特征，這表明殘差分布并不完全服從正態(tài)分布。采用Shapiro等[28]提出的W′正態(tài)檢驗(yàn)方法，計(jì)算了統(tǒng)計(jì)量W′=0.991，達(dá)到了統(tǒng)計(jì)顯著性95%信度水平，認(rèn)為殘差通過了正態(tài)性檢驗(yàn)，服從正態(tài)分布。Montgomery[29]指出使用GAM時(shí)，樣本低于100時(shí)有可能出現(xiàn)類似圖3的現(xiàn)象。

圖3 模型殘差正態(tài)概率圖和殘差直方圖Fig.3 Normal probability plot and histogram of model residuals

為進(jìn)一步驗(yàn)證模型的適用性，探討了模型殘差散點(diǎn)分布以及實(shí)測(cè)流量與模擬流量散點(diǎn)的擬合情況。由圖4(a)可知，殘差大部分處于-50～ 50之間，分布基本合理。從圖4(b)可知，其值大致呈現(xiàn)1∶1線性分布的趨勢(shì)，流量存在一個(gè)異常大值(>500)也得到了較為準(zhǔn)確的模擬，表明模型的擬合效果較好。

圖4 殘差散點(diǎn)圖及流量的實(shí)測(cè)值和模擬值散點(diǎn)分布Fig.4 Residual scatter plot and scatter plot of observed runoff versus predicted runoff

綜上所述，可以認(rèn)定GAM用于本文的徑流模擬是合適的。

2.3 模型預(yù)測(cè)效果評(píng)估

采用留一法交叉驗(yàn)證對(duì)GAM

和GLM的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比分析，其中樣本數(shù)據(jù)為1952—2014年。根據(jù)留一法交叉驗(yàn)證的原理，每次均選取其中的62 a作為訓(xùn)練集，剩下的1 a作為測(cè)試集驗(yàn)證模型的優(yōu)劣，循環(huán)往復(fù)，共進(jìn)行63次模擬與預(yù)測(cè)。

圖5為實(shí)測(cè)值與2種模型預(yù)測(cè)值的對(duì)比，其中實(shí)測(cè)值與GLM預(yù)測(cè)序列的相關(guān)性為0.267，通過了信度95%的顯著性水平；實(shí)測(cè)值與GAM預(yù)測(cè)序列的相關(guān)系數(shù)為0.407，通過了信度99.99%的顯著性檢驗(yàn)。因此，相較于GLM，GAM的預(yù)測(cè)能力也有了顯著提升。在一些流量異常年份的預(yù)測(cè)上，如1957年、1968年、2001年、2008年幾個(gè)流量偏多年份，GAM的預(yù)測(cè)效果相較于GLM均有明顯提升，而在1985年和1998年2個(gè)GLM預(yù)測(cè)值誤差顯著偏大的年份，GAM的預(yù)測(cè)結(jié)果也有了很大的改善。雖然也有部分年份GLM預(yù)測(cè)更接近實(shí)測(cè)值，但從整體趨勢(shì)及一些流量異常的年份看，GAM的預(yù)測(cè)能力均高于GLM。

圖5 2種模型的預(yù)測(cè)值與實(shí)測(cè)值對(duì)比Fig.5 Comparison between prediction values of GAM and GLM and observed values

從圖6可以看到在10%的相對(duì)誤差絕對(duì)值范圍內(nèi)，GAM預(yù)報(bào)正確年份為13 a，與GLM相比少1 a，而在10%～20%以及20%～30%的誤差范圍內(nèi)GAM預(yù)報(bào)準(zhǔn)確的年份比GLM分別多3 a和5 a。而當(dāng)相對(duì)誤差絕對(duì)值大于40%時(shí)，GLM預(yù)報(bào)的年份要明顯多于GAM。因此，在頻數(shù)分布上也可以看出，GAM模型的預(yù)報(bào)效果明顯優(yōu)于GLM的預(yù)報(bào)效果。

圖6 GAM和 GLM預(yù)測(cè)值相對(duì)誤差頻數(shù)分布Fig.6 Frequency distribution of prediction relative error of GAM and GLM

如果以預(yù)測(cè)相對(duì)誤差絕對(duì)值小于等于30%為合格的預(yù)測(cè)結(jié)果，GAM和GLM的預(yù)測(cè)合格率分別為85.5%和79.0%。表明GAM的預(yù)測(cè)合格率高于GLM，預(yù)測(cè)效果更好。

文獻(xiàn)[15]使用了5種模型(時(shí)間序列方法、周期均值疊加法、逐步回歸、嶺回歸和主成分回歸)對(duì)2004年和2005年分別進(jìn)行了預(yù)測(cè)分析，從6—8月平均相對(duì)誤差看，這2 a平均相對(duì)誤差分別是18.0%和142.0%，而GAM對(duì)這兩年預(yù)測(cè)的相對(duì)誤差分別是9.0%和68.8%，預(yù)測(cè)精度均明顯提高。這也印證了考慮前期海溫因子的非線性作用后，模型模擬及預(yù)測(cè)能力均有明顯提高。因此，GAM能夠?qū)鹘y(tǒng)線性回歸模型的預(yù)報(bào)效果加以改進(jìn)。

3 結(jié) 語

在前期海溫因子中，挑選了4個(gè)非線性因子和9個(gè)線性因子，構(gòu)建的GAM對(duì)洪家渡的夏季徑流預(yù)測(cè)非常有效，通過AIC、RMSE、MAE、LEPS、相關(guān)性5種評(píng)價(jià)指標(biāo)評(píng)估，GAM都優(yōu)于GLM。

對(duì)洪家渡夏季徑流的預(yù)測(cè)評(píng)估結(jié)果表明， GAM的預(yù)測(cè)效果比GLM明顯提高，從相對(duì)誤差絕對(duì)值的頻數(shù)、極值預(yù)報(bào)誤差和預(yù)測(cè)合格率等3個(gè)方面的分析表明，GAM的誤差范圍更小，極值預(yù)測(cè)的精度和預(yù)測(cè)合格率均比較高。

徑流的中長期預(yù)測(cè)對(duì)水庫調(diào)度和水資源管理有重要作用，研究高精度的徑流預(yù)測(cè)模型意義重大。本文的嘗試盡管提高了洪家渡夏季徑流的預(yù)測(cè)能力，但是，今后通過因子影響機(jī)理分析和消除因子之間的共線性等問題，選取更好的預(yù)測(cè)指標(biāo)，或許可以進(jìn)一步提高GAM的模擬精度。