程忠良，劉 勇，高 成，胡 健

(1. 河海大學(xué)水文水資源學(xué)院，南京 210098；2. 南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210029；3. 中國(guó)科學(xué)院 陸地水循環(huán)及地表過(guò)程重點(diǎn)實(shí)驗(yàn)室，北京 100101)

0 引 言

徑流預(yù)報(bào)是水文學(xué)的重要應(yīng)用領(lǐng)域之一，也是水庫(kù)運(yùn)行調(diào)度和防汛、抗旱以及水資源應(yīng)急調(diào)度等的前提[1]。長(zhǎng)期徑流預(yù)報(bào)因其預(yù)見期較長(zhǎng)、預(yù)報(bào)精度不高，難以滿足人們對(duì)社會(huì)經(jīng)濟(jì)生產(chǎn)安排的需求，一直以來(lái)都是水文領(lǐng)域研究的難點(diǎn)之一。水文預(yù)測(cè)方法眾多，通常大致可分為過(guò)程驅(qū)動(dòng)和數(shù)據(jù)驅(qū)動(dòng)兩大類[2]。其中過(guò)程驅(qū)動(dòng)法是基于產(chǎn)匯流機(jī)制建立的模型，是徑流預(yù)測(cè)的一個(gè)發(fā)展方向。但由于徑流受到氣候氣象、下墊面、人類活動(dòng)等多方面不確定性因素的影響[3]，其形成機(jī)理和規(guī)律尚未被完全掌握，預(yù)報(bào)精度不高，使得該方法的應(yīng)用十分困難。隨著數(shù)據(jù)獲取能力和計(jì)算能力的發(fā)展，數(shù)據(jù)驅(qū)動(dòng)模型在水文預(yù)報(bào)中的應(yīng)用也愈加廣泛，神經(jīng)網(wǎng)絡(luò)[4]、支持向量機(jī)[5]等預(yù)測(cè)方法在實(shí)際應(yīng)用中也取得了一定的成果。近年來(lái)，通過(guò)發(fā)掘未來(lái)徑流與前期降雨、海表溫度、大氣環(huán)流指數(shù)等大尺度氣候-水文變量間的關(guān)聯(lián)關(guān)系，建立基于物理成因背景的徑流預(yù)報(bào)模型是目前研究的重要方向[6-9]。馮小沖[6]、劉勇[7]等通過(guò)對(duì)影響丹江口水庫(kù)入庫(kù)徑流的氣候和水文氣象物理因子的分析，篩選出與徑流密切相關(guān)的物理因子作為預(yù)報(bào)因子，分別采用逐步回歸和神經(jīng)網(wǎng)絡(luò)的方法構(gòu)建月徑流預(yù)報(bào)模型，取得了一定的成果。在氣候變化和高強(qiáng)度人類活動(dòng)的背景下，全球海洋和大氣環(huán)流均發(fā)生了顯著的年代際變化[10]，水庫(kù)來(lái)水量的變化也隨之加劇，單獨(dú)采用定性或定量的回歸預(yù)報(bào)模型存在預(yù)報(bào)誤差偏大和精度偏低的問(wèn)題。為進(jìn)一步提高預(yù)報(bào)的精度和可靠性，本文擬在前人研究的基礎(chǔ)上，采用定性定量相結(jié)合的預(yù)報(bào)方法，從影響丹江口水庫(kù)長(zhǎng)期徑流的物理成因出發(fā)選擇預(yù)報(bào)因子，采用AIC準(zhǔn)則篩選關(guān)鍵因子，劃分徑流級(jí)別，利用馬氏距離判別分析的分類判別優(yōu)勢(shì)，構(gòu)建長(zhǎng)期徑流分級(jí)預(yù)報(bào)模型，并與逐步回歸模型作比較。

1 研究思路及方法

1.1 研究思路

首先從物理成因出發(fā)，考慮影響丹江口水庫(kù)入庫(kù)徑流的海溫、環(huán)流等因素，采用單相關(guān)系數(shù)法，計(jì)算預(yù)報(bào)因子與預(yù)報(bào)徑流之間的相關(guān)系數(shù)r，并作0.05信度水平下的檢驗(yàn)，初選出基礎(chǔ)預(yù)報(bào)因子集。然后，將資料序列劃分為模擬期和檢驗(yàn)期，模擬期的徑流按照距平要素值劃分為豐、平、枯三個(gè)級(jí)別，采用AIC準(zhǔn)則進(jìn)一步篩選出關(guān)鍵因子集，構(gòu)建分段多元線性回歸方程，對(duì)預(yù)報(bào)對(duì)象進(jìn)行模擬；檢驗(yàn)期時(shí)，采用馬氏距離判別法根據(jù)預(yù)報(bào)因子對(duì)預(yù)報(bào)徑流先做出定性判別，再依據(jù)定性判別類別代入相應(yīng)的線性回歸方程，做出定量預(yù)報(bào)，并將該預(yù)報(bào)模型與逐步回歸模型的模擬和預(yù)報(bào)結(jié)果進(jìn)行對(duì)比，從預(yù)報(bào)精度和穩(wěn)定性兩個(gè)方面去比較兩種模型。

1.2 研究方法

1.2.1 單相關(guān)系數(shù)法

在水文中長(zhǎng)期預(yù)測(cè)中常采用相關(guān)系數(shù)來(lái)考察預(yù)報(bào)因子與預(yù)報(bào)對(duì)象之間是否線性相關(guān)，并以此作為因子挑選的依據(jù)[3]。單相關(guān)系數(shù)的公式為：

(1)

1.2.2 AIC準(zhǔn)則

因子個(gè)數(shù)p的選取對(duì)模型的預(yù)報(bào)精度和穩(wěn)定性具有很大的影響。取較小的p，擬合程度較差；取較大的p，則受偶然變化影響太大，甚至出現(xiàn)過(guò)度擬合的現(xiàn)象，預(yù)報(bào)效果會(huì)受到影響。AIC是衡量統(tǒng)計(jì)模型擬合優(yōu)良性的一種標(biāo)準(zhǔn)，由日本統(tǒng)計(jì)學(xué)家赤池弘次在1974年提出，它建立在熵的概念上，提供了權(quán)衡估計(jì)模型復(fù)雜度和擬合數(shù)據(jù)優(yōu)良性的標(biāo)準(zhǔn)?，F(xiàn)定義一個(gè)準(zhǔn)則函數(shù)為：

(2)

式中：σ2(p)為誤差的方差；p為因子個(gè)數(shù)；n為樣本個(gè)數(shù)。

其中第一項(xiàng)代表的是擬合優(yōu)度，第二項(xiàng)代表了增加因子后的懲罰，權(quán)衡兩者，選取最小的AIC(p)作為合理的因子個(gè)數(shù)p。利用AIC準(zhǔn)則篩選出關(guān)鍵預(yù)報(bào)因子集，用于建立回歸方程。

1.2.3 馬氏距離判別分析

馬氏距離是由印度統(tǒng)計(jì)學(xué)家馬哈拉諾比斯(Mahalanobis)提出的，表示數(shù)據(jù)的協(xié)方差距離。馬氏距離判別分析法是根據(jù)觀測(cè)到的樣本的若干數(shù)量特征對(duì)新獲得的樣本進(jìn)行歸類、識(shí)別，判別其所屬類型的一種有效計(jì)算樣本集相似度的統(tǒng)計(jì)分析方法。相比于歐氏距離判別方法，馬氏距離的優(yōu)勢(shì)在于不受量綱的影響，兩點(diǎn)間的馬氏距離與原始數(shù)據(jù)的測(cè)量無(wú)關(guān)，可排除變量間相關(guān)性的干擾。該方法的主要思想是比較樣本到各個(gè)總體的馬氏距離，然后將其判給馬氏距離最近的那個(gè)總體。

設(shè)有k個(gè)m維總體G1，G2，…，Gk，均值向量分別為μ1，μ2，…，μk，協(xié)方差矩陣分別為∑1，∑2，…，∑k，則樣本X到各組的馬氏距離為：

D2(X,Gα)=(X-μα)T∑α-1(X-μα)

(3)

α=1，2，…，k

判別規(guī)則為：若D2(X，Gi)=min1≤α≤kD2(X，Gα)則X∈Gi。

針對(duì)實(shí)際問(wèn)題，當(dāng)μ1，μ2，…，μk和∑1，∑2，…，∑k均未知時(shí)，可以通過(guò)相應(yīng)的樣本估計(jì)值來(lái)替代，馬氏距離的具體計(jì)算步驟如下。

(4)

(5)

(6)

采用馬氏距離判別法，計(jì)算樣本到各個(gè)總體的馬氏距離，利用馬氏距離最小準(zhǔn)則對(duì)預(yù)報(bào)期徑流的豐、平、枯類別做出定性預(yù)報(bào)。

1.2.4 多元線性回歸

多元線性回歸是用于研究一個(gè)隨機(jī)變量和多個(gè)變量之間相關(guān)關(guān)系的方法，利用多個(gè)因子X1，X2，…，Xm與對(duì)象Y建立多元回歸方程：

Y=b0+b1X1+b2X2+…+bmXm

(7)

其中b0，b1，…，bm為回歸系數(shù)，多元線性回歸方程的回歸系數(shù)采用最小二乘法來(lái)確定。根據(jù)預(yù)報(bào)因子和預(yù)報(bào)對(duì)象之間的關(guān)系可以建立相應(yīng)的回歸方程，做出定量預(yù)報(bào)。

1.3 模型精度評(píng)定方法

模型的精度和穩(wěn)定性采用《水文情報(bào)預(yù)報(bào)規(guī)范》(GB/T22482-2008)[11]中的確定性系數(shù)和平均絕對(duì)誤差來(lái)評(píng)定，確定性系數(shù)DC和平均絕對(duì)誤差e的計(jì)算公式見式(8)和式(9)。

(8)

(9)

2 丹江口9月入庫(kù)徑流預(yù)報(bào)

2.1 研究區(qū)域概況

丹江口水庫(kù)位于漢江中上游，丹江與漢江干流匯合以下800 m處，其地理位置處于東經(jīng)106°12′～111°26′，北緯31°24′～34°11′之間，是南水北調(diào)中線工程的水源地。漢江丹江口水庫(kù)具有防洪、發(fā)電、灌溉、航運(yùn)和養(yǎng)殖等五大功能，水庫(kù)多年平均入庫(kù)水量為394.8 億m3，控制流域面積95 217 km2,是我國(guó)大型綜合利用水庫(kù)之一。

2.2 預(yù)報(bào)因子初選

利用收集到的丹江口水庫(kù)1952-2008年的9月入庫(kù)徑流資料以及北太平洋逐月平均海溫、北半球100 hPa和500 hPa逐月平均高度場(chǎng)以及74項(xiàng)逐月環(huán)流特征量資料，分別對(duì)各要素場(chǎng)的因子進(jìn)行相關(guān)性普查，并通過(guò)0.05的顯著性檢驗(yàn)。初步篩選出北半球100 hPa逐月平均高度場(chǎng)因子5項(xiàng)、北半球500 hPa逐月平均高度場(chǎng)因子7項(xiàng)、北太平洋海溫場(chǎng)因子5項(xiàng)、74項(xiàng)環(huán)流特征量因子6項(xiàng)，建立徑流預(yù)報(bào)基礎(chǔ)因子集。

2.3 基于馬氏距離判別分析的長(zhǎng)期徑流分級(jí)預(yù)報(bào)

把丹江口水庫(kù)1952-2008年9月入庫(kù)徑流量資料分為模擬期(1952-2000年)和檢驗(yàn)期(2001-2008年)，并根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》(GB/T22482-2008)[11]，將丹江口水庫(kù)模擬期的9月入庫(kù)徑流按照距平值劃分為豐、平、枯三段，詳細(xì)劃分標(biāo)準(zhǔn)見表1。

表1 豐、平、枯級(jí)別劃分標(biāo)準(zhǔn)Tab.1 standard of division of high, flat and dry

采用AIC準(zhǔn)則進(jìn)行關(guān)鍵預(yù)報(bào)因子的遴選。篩選出的關(guān)鍵因子集見表2。

從篩選的因子集來(lái)看，預(yù)報(bào)因子對(duì)于不同徑流級(jí)別的預(yù)報(bào)影響存在差異，海溫因子對(duì)偏豐徑流的預(yù)報(bào)影響較大，而100 hPa和500 hPa位勢(shì)高度對(duì)偏枯徑流的預(yù)報(bào)影響較大，未發(fā)現(xiàn)對(duì)豐、平、枯均有良好預(yù)報(bào)影響的因子。

表2 AIC準(zhǔn)則篩選的關(guān)鍵因子集Tab.2 key factor sets selected by AIC

根據(jù)AIC準(zhǔn)則篩選出的關(guān)鍵影響因子，利用多元線性回歸法，建立的分段線性回歸方程如下：

(10)

根據(jù)建立的回歸方程，可對(duì)模擬期1952-2000年的9月入庫(kù)徑流進(jìn)行模擬預(yù)報(bào)。

根據(jù)篩選出的關(guān)鍵因子集，利用馬氏距離判別法進(jìn)行判別分析，判別分析結(jié)果見表3。

表3 馬氏距離判別計(jì)算表Tab.3 Mahalanobis distance discriminate calculation table

訓(xùn)練期1952-2000年的49年中有43年判定合格，合格率為87.8%；檢驗(yàn)期2001-2008年的8年中有7年判定合格，合格率為87.5%，具有較高精度。對(duì)于檢驗(yàn)期2001-2008年，根據(jù)馬氏距離判別分析的類別，代入相應(yīng)的回歸方程，做出預(yù)報(bào)。

2.4 逐步回歸法預(yù)報(bào)模型

根據(jù)預(yù)報(bào)初選因子集，采用逐步回歸分析的方法進(jìn)一步篩選和確定關(guān)鍵影響因子集，篩選出的關(guān)鍵影響因子集見表4。

利用篩選出的關(guān)鍵因子集構(gòu)建回歸方程，預(yù)報(bào)方程見式(11)。

Y=99.96+0.206X1-0.180X2-0.142X3-0.397X4+

20.1X5+3.60X6+0.799X7+ 5.03X8-6.70X9-2.23X10

(11)

3 預(yù)報(bào)成果對(duì)比分析

利用基于馬氏距離判別的長(zhǎng)期徑流預(yù)報(bào)模型(方法1)和逐步回歸預(yù)報(bào)模型(方法2)分別對(duì)丹江口水庫(kù)9月入庫(kù)徑流量進(jìn)行模擬預(yù)報(bào)，模擬期預(yù)報(bào)成果見圖1。

表4 逐步回歸分析篩選的關(guān)鍵因子集Tab.4 key factor sets selected by stepwise regression analysis

圖1 丹江口水庫(kù)9月入庫(kù)徑流量實(shí)測(cè)值與模擬值Fig.1 The measured value and simulated value of the Danjiangkou reservoir in September

從模擬結(jié)果來(lái)看，根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》(GB/T22482-2008)[11]中評(píng)定中長(zhǎng)期預(yù)報(bào)精度的方案作為預(yù)報(bào)精度模型預(yù)報(bào)的評(píng)價(jià)標(biāo)準(zhǔn)：對(duì)于定量預(yù)報(bào)，水位(流量)按多年變幅的10%、其他要素按20%，要素極值的出現(xiàn)時(shí)間按多年變幅的30%作為許可誤差。采用多年變幅的10%(22.73 億m3)作為許可誤差，方法1的49年模擬期中有46年滿足要求，合格率為93.8%，方法2的49年模擬期中有44年滿足要求，合格率為89.8%。兩種方法的確定性系數(shù)DC和平均絕對(duì)誤差e見表5。

表5 模型確定性系數(shù)和平均絕對(duì)誤差統(tǒng)計(jì)表Tab.5 statistical table of deterministic coefficientand mean absolute error

在檢驗(yàn)期2001-2008年的8年中，兩種方法的實(shí)測(cè)值與預(yù)報(bào)值的誤差見表6。

根據(jù)多年變幅的10%作為許可誤差，方法1中只有2003年的預(yù)報(bào)誤差超過(guò)允許值，預(yù)報(bào)的合格率為87.5%；方法2中的2003年和2007年的誤差超過(guò)允許值，合格率為75%。

表6 檢驗(yàn)期實(shí)測(cè)值與預(yù)報(bào)值誤差統(tǒng)計(jì)表Tab.6 statistical table of measured and predictedvalues in the test period

4 結(jié) 語(yǔ)

(1)預(yù)報(bào)結(jié)果表明，基于馬氏距離判別的分級(jí)預(yù)報(bào)模型采用定性定量預(yù)報(bào)相結(jié)合的方法，預(yù)報(bào)的精度和穩(wěn)定性優(yōu)于采用單一定量預(yù)報(bào)的逐步回歸模型。

(2)利用AIC準(zhǔn)則通過(guò)綜合考慮模型精度和復(fù)雜度來(lái)確定因子個(gè)數(shù)，對(duì)于提高預(yù)報(bào)精度和穩(wěn)定性有一定的幫助，提高預(yù)報(bào)精度關(guān)鍵在于提高因子的代表性而不僅僅是增加因子的數(shù)量。

(3)基于物理成因出發(fā)篩選出的預(yù)報(bào)因子對(duì)提高預(yù)報(bào)精度有較好的效果，今后可加強(qiáng)預(yù)報(bào)因子之間的關(guān)聯(lián)性分析以及預(yù)報(bào)因子與預(yù)報(bào)對(duì)象之間的成因機(jī)理分析研究，進(jìn)一步提高預(yù)報(bào)因子的相關(guān)性和代表性，從而減小預(yù)報(bào)的誤差。