基于混合回歸模型的沿海平原河網(wǎng)洪水位預(yù)報(bào)研究

周 煥，閃麗潔，揭夢(mèng)璇，湯 斌

（浙江省水利水電勘測(cè)設(shè)計(jì)院，浙江 杭州 310002）

1 研究背景

浙江省位于我國(guó)東南沿海、太平洋西岸的亞熱帶季風(fēng)區(qū)，是一個(gè)洪澇臺(tái)災(zāi)害易發(fā)多發(fā)的地區(qū)。東部沿海平原既是浙江省經(jīng)濟(jì)發(fā)展主導(dǎo)地區(qū)，也是人口和生產(chǎn)力要素聚集區(qū)，多處在河流中下游，背山濱海，上承山洪，下受潮汐頂托，地勢(shì)低平，排水不暢。近年來(lái)，極端天氣多發(fā)，風(fēng)、暴、潮不利組合幾率增加，競(jìng)賽式圩區(qū)建設(shè)導(dǎo)致外江洪水位不斷抬升，加上部分區(qū)域地面沉降、入海排水口門外延等諸多因素，造成平原排澇形勢(shì)“雪上加霜”，對(duì)平原感潮河網(wǎng)城市防洪排澇工作帶來(lái)巨大的挑戰(zhàn)［1-3］。

平原感潮河網(wǎng)城市防洪排澇的關(guān)鍵在于控制市區(qū)境內(nèi)主要河道的水位。相比于山區(qū)的單一或樹(shù)狀河道，平原河網(wǎng)水文條件錯(cuò)綜復(fù)雜，由此帶來(lái)預(yù)報(bào)模型的方程組離散和求解更加困難，這是多年來(lái)人們研究河網(wǎng)問(wèn)題的難點(diǎn)［4］。以往研究中，主要采用水動(dòng)力模型，涉及控制方程組簡(jiǎn)化、方程組離散和求解、初邊條件確定等一系列問(wèn)題［5-7］。此類方法需要詳細(xì)的河道地形資料及特征參數(shù)，且計(jì)算收斂穩(wěn)定性受輸入邊界影響較大，規(guī)模龐大、時(shí)間長(zhǎng)。近年來(lái)隨著物聯(lián)網(wǎng)技術(shù)的發(fā)展和數(shù)據(jù)感知技術(shù)的成熟，多元逐步回歸模型、人工神經(jīng)網(wǎng)絡(luò)模型等基于大數(shù)據(jù)的數(shù)理統(tǒng)計(jì)方法快速發(fā)展，在氣象、水文等方面得到了大量的應(yīng)用［8-11］。都宏博［12］采用門限回歸和BP人工神經(jīng)網(wǎng)絡(luò)模型建立了感潮河段水位預(yù)報(bào)模型，結(jié)果表明構(gòu)建模型比較符合感潮河段水位變化的物理規(guī)律，適用于長(zhǎng)江流域感潮河段水位預(yù)報(bào)。王竹等［13］將智能算法與半分布式水文模型相結(jié)合，構(gòu)建了一種半分布式BP神經(jīng)網(wǎng)絡(luò)模型，并應(yīng)用在遼寧省大伙房水庫(kù)洪水預(yù)報(bào)工作中，取得了理想的預(yù)報(bào)效果。然而數(shù)理統(tǒng)計(jì)方法在平原河網(wǎng)洪水預(yù)報(bào)工作中應(yīng)用較少。

本文以浙江沿海溫瑞平原的9個(gè)特征斷面為研究對(duì)象，以各斷面2005年5號(hào)“海棠”臺(tái)風(fēng)、2007年13號(hào)“韋帕”臺(tái)風(fēng)、2009年8號(hào)“莫拉克”臺(tái)風(fēng)、2013年23號(hào)“菲特”臺(tái)風(fēng)等6場(chǎng)洪水過(guò)程為依據(jù)，嘗試采用在水文產(chǎn)匯流預(yù)報(bào)模擬成果的基礎(chǔ)上利用數(shù)理統(tǒng)計(jì)方法解決問(wèn)題的思路，利用產(chǎn)匯流模型計(jì)算求得各水文分區(qū)洪水過(guò)程，與下游潮位過(guò)程一起形成預(yù)報(bào)因子庫(kù)，并在此基礎(chǔ)上針對(duì)各特征斷面構(gòu)建相應(yīng)的混合回歸模型，從而進(jìn)行溫瑞平原洪水位預(yù)報(bào)研究。

2 研究流域和數(shù)據(jù)資料

溫瑞平原是浙江省東南沿海甌江和飛云江之間多塊平原的總稱。溫瑞平原總集水面積948 km2。地形總體上西部及西南部為山區(qū)、東部為平原，北臨甌江、南界飛云江、東瀕東海，地理位置如圖1所示。區(qū)域整體地形呈西高東低，山區(qū)無(wú)調(diào)蓄水庫(kù)、山腳無(wú)集水溝，上游洪水直接排入城區(qū)或低地，城區(qū)易成為上游周邊山區(qū)和較高地勢(shì)來(lái)水的蓄水區(qū)［14］。平原河網(wǎng)水系分布較密，其中主要水系為溫瑞塘河，北起溫州，南至瑞安市城區(qū)，是溝通甌江、飛云江兩大水系的主要內(nèi)河［15］。

圖1 溫瑞平原流域分布圖

經(jīng)分析得知，溫瑞平原排澇格局為向北排甌江、向南排飛云江、向東排東海。為了更加準(zhǔn)確地模擬特征斷面水位，本文根據(jù)平原河網(wǎng)內(nèi)水的排流情況，將溫瑞平原流域劃分為3個(gè)片區(qū)，分別為鹿城甌海片、龍灣片和瑞安片，并對(duì)不同片區(qū)內(nèi)的斷面水位分別進(jìn)行模擬。所選擇的率定期為1999年9月3日—6日、2005年7月18日—21日、2007年9月18日—21日和2009年8月7日—10日共4場(chǎng)洪水，檢驗(yàn)期為2010年7月24日—27日和2013年10月6日—9日共2場(chǎng)洪水。所用的資料包括6場(chǎng)洪水對(duì)應(yīng)的流域內(nèi)雨量站降雨資料、9個(gè)特征斷面水位資料、臨海閘門閘前水位資料、閘門參數(shù)資料以及瑞安站潮位資料等，時(shí)間尺度為小時(shí)。

3 研究方法

與山區(qū)型河道匯流速度快、匯流方向單一的特點(diǎn)不同，溫瑞平原地勢(shì)平坦，水系四通八達(dá)，特征斷面洪水過(guò)程很大程度受周圍片區(qū)水位變化的影響，因此平原地區(qū)洪水位預(yù)報(bào)更加復(fù)雜。參考傳統(tǒng)的水動(dòng)力模型如MIKE11的建模思路，首先對(duì)整個(gè)溫瑞水系區(qū)域河道現(xiàn)狀水力條件進(jìn)行分析，結(jié)合水文站和雨量站分布劃分為子流域，建立產(chǎn)匯流模型模擬各子流域的降雨徑流過(guò)程，作為溫瑞平原混合回歸模型的預(yù)報(bào)因子庫(kù)，并在此基礎(chǔ)上構(gòu)建混合回歸模型，即首先分別構(gòu)建門限回歸模型和多元逐步回歸模型，再采用多元線性回歸模型進(jìn)行集合預(yù)報(bào)，從而實(shí)現(xiàn)對(duì)溫瑞平原洪水預(yù)報(bào)工作。

3.1 水文分區(qū)產(chǎn)匯流模型構(gòu)建 根據(jù)流域特性將研究區(qū)域劃分多個(gè)水文分區(qū)，根據(jù)站點(diǎn)觀測(cè)的雨量計(jì)算各水文分區(qū)的流量過(guò)程，為混合回歸綜合模型提供預(yù)報(bào)因子庫(kù)。依據(jù)河流水系、分水嶺及下墊面情況等要素，本研究區(qū)域可細(xì)化39個(gè)山區(qū)分區(qū)和91個(gè)平原分區(qū)。采用泰森多邊形法由8個(gè)雨量站點(diǎn)的雨量加權(quán)計(jì)算得到各分區(qū)的面雨量，再分別對(duì)山區(qū)和平原分區(qū)采用不同的模型計(jì)算產(chǎn)匯流過(guò)程。各分區(qū)邊界范圍如圖2所示。

圖2 研究區(qū)域水文分區(qū)圖

山區(qū)分區(qū)采用蓄滿產(chǎn)流模型進(jìn)行產(chǎn)流計(jì)算，后根據(jù)分區(qū)集水面積大小，選用不同的計(jì)算方法，集水面積大于50 km2的采用瞬時(shí)單位線法，小于50 km2的采用“浙江省推理公式法”，具體計(jì)算方法參考論文［16］。平原分區(qū)按不同地類分別推算產(chǎn)水過(guò)程：水面按水量平衡方程由降雨扣除水面蒸發(fā)推求產(chǎn)水量；水田由降雨扣除水稻蒸騰及水田下滲并考慮水田最大持水深度推求產(chǎn)水量；旱地根據(jù)前期土壤含水量由降雨扣除一定的損失及下滲，并考慮旱地最大持水深度從而推求產(chǎn)水量；其他建成區(qū)則采用徑流系數(shù)法由降雨推求產(chǎn)流過(guò)程。

3.2 預(yù)報(bào)因子的確定 預(yù)報(bào)因子篩選準(zhǔn)確、全面與否，是對(duì)準(zhǔn)確模擬洪水位工作的前提［17］。結(jié)合流域地理、氣候、洪水特征等特點(diǎn)，采取物理成因分析與數(shù)理統(tǒng)計(jì)分析相結(jié)合的思路，從對(duì)洪水過(guò)程的貢獻(xiàn)率、影響程度及其獲取難易程度等方面綜合考慮，識(shí)別洪水位的主要影響因子，為模型構(gòu)建奠定必要的基礎(chǔ)［18-19］。

為了得到更為合理準(zhǔn)確的預(yù)報(bào)因子組合，本文采用相關(guān)系數(shù)與逐步回歸耦合的方法分兩步進(jìn)行篩選工作［20］。首先選取足夠數(shù)量的影響因子作為預(yù)報(bào)因子庫(kù)，在成因分析的基礎(chǔ)上分析洪水來(lái)源，同時(shí)分析特征斷面周圍的其他影響因素。其次，采用相關(guān)系數(shù)法從預(yù)報(bào)因子庫(kù)中篩選出與斷面洪水位過(guò)程相關(guān)性較強(qiáng)的因子，作為待選預(yù)報(bào)因子組合。最后，通過(guò)逐步回歸模型對(duì)待選因子進(jìn)一步篩選，得到最終的預(yù)報(bào)因子組合，作為混合回歸模型的輸入項(xiàng)進(jìn)行洪水位預(yù)報(bào)。

3.3 混合回歸模型構(gòu)建 混合回歸模型是在模擬洪水變化過(guò)程尋求相關(guān)鏈接因子的一種自適應(yīng)統(tǒng)計(jì)學(xué)技術(shù)，主要由門限回歸模型和多元逐步回歸模型組成，思路為分別構(gòu)建門限回歸模型和多元逐步回歸模型，并在此基礎(chǔ)上采用多元線性回歸法建立集合預(yù)報(bào)模型。

3.3.1 門限回歸模型 門限回歸模型由英籍華人H.Tong（湯家豪）博士首創(chuàng)提出，將傳統(tǒng)的全局線性逼近分成幾段進(jìn)行線性回歸逼近，分割由“門限”變量來(lái)控制［21-22］。其實(shí)質(zhì)上是把預(yù)報(bào)問(wèn)題按狀態(tài)空間的取值進(jìn)行分類，用分段的線性回歸模式來(lái)描述總體非線性預(yù)報(bào)問(wèn)題［23］。本文采用二分割門限回歸模型，將待模擬洪水位序列及其它預(yù)報(bào)因子序列按最優(yōu)分割點(diǎn)分為兩段，分別進(jìn)行多元回歸計(jì)算。門限回歸模型可表示為：

3.3.2 多元逐步回歸模型 多元逐步回歸模型是在多元線性回歸分析的基礎(chǔ)上派生的一種建立最優(yōu)多元線性回歸方程的算法技巧，核心是利用求解求逆緊湊變化法和雙檢驗(yàn)法，從大量因子中選擇重要有效的因子建立最優(yōu)回歸方程［24］。目前，我國(guó)廣泛應(yīng)用的逐步回歸算法，其實(shí)質(zhì)是用高斯-亞當(dāng)消去法解正規(guī)方程式，巧妙地把選取重要因子的思想結(jié)合進(jìn)去，大大提高了計(jì)算效率［17］。

3.4 評(píng)價(jià)指標(biāo) 本文所采用的洪水預(yù)報(bào)精度評(píng)定項(xiàng)目包括洪峰水位、洪峰水位出現(xiàn)時(shí)間和洪水過(guò)程。選擇確定性系數(shù)Dc和合格率QR為洪水預(yù)報(bào)精度評(píng)價(jià)指標(biāo)。依據(jù)《水文情報(bào)設(shè)計(jì)規(guī)范》中定義，計(jì)算公式如下：

4 結(jié)果與分析

4.1 預(yù)報(bào)因子篩選 首先，在成因分析的基礎(chǔ)上分析洪水來(lái)源。溫瑞平原洪水受上游山區(qū)洪水下泄、區(qū)間暴雨和下游潮位頂托的共同影響。其中上游山區(qū)洪水下泄、區(qū)間暴雨均可通過(guò)產(chǎn)匯流模型計(jì)算各分區(qū)產(chǎn)匯流過(guò)程進(jìn)行體現(xiàn)。另外，溫瑞平原沿海地區(qū)分布較多擋潮閘和排（退）水閘，為了在排洪水的同時(shí)防止潮水倒灌，閘門會(huì)隨著外海潮位的漲落而不斷啟閉，閘上水位亦會(huì)隨之發(fā)生變化。通過(guò)分析發(fā)現(xiàn)，康華僑東、西山、水心殿、狀元等斷面洪水過(guò)程均會(huì)不同程度受到沿海閘門啟閉的影響，因此斷面附近閘門閘上水位亦可作為預(yù)報(bào)因子，與各水文分區(qū)產(chǎn)匯流過(guò)程共同形成預(yù)報(bào)因子庫(kù)。其次，通過(guò)相關(guān)系數(shù)法逐時(shí)計(jì)算因子與洪水位序列的相關(guān)系數(shù)，選出相關(guān)系數(shù)絕對(duì)值最大的10～15個(gè)因子為待選預(yù)報(bào)因子。最后，通過(guò)逐步回歸模型對(duì)待選因子進(jìn)一步篩選，得到最終的預(yù)報(bào)因子組合，作為混合回歸模型的輸入項(xiàng)進(jìn)行洪水位預(yù)報(bào)。

需要注意的是，閘上水位過(guò)程需要通過(guò)堰流公式計(jì)算得出。首先根據(jù)流域?qū)嶋H情況，匯總水文邊界入分區(qū)產(chǎn)流過(guò)程，對(duì)區(qū)域總流量過(guò)程進(jìn)行模擬，并根據(jù)各閘門大小以及所處位置，模擬出總流量過(guò)程對(duì)各閘口斷面流量的分配關(guān)系，從而進(jìn)一步得到各閘口斷面的流量過(guò)程。模型將閘門處至閘上斷面間的河道概化為一個(gè)水庫(kù)，依據(jù)水量平衡原理，采用堰流流量公式計(jì)算下泄水量和閘上水位過(guò)程。為了更加精確地模擬閘前水位和斷面入流過(guò)程的變化，將時(shí)間尺度由1 h縮短至15 s，對(duì)應(yīng)的流量過(guò)程和下邊界潮位過(guò)程由插值所得。

結(jié)合溫瑞平原各斷面的地理位置、洪水來(lái)源等因素，最終篩選出的預(yù)報(bào)因子組合如表1所示。可以看出，西山、康華僑東和水心殿斷面與勤奮水閘較近，斷面水位過(guò)程與勤奮閘的閘前水位過(guò)程存在較大的相關(guān)性；永強(qiáng)和沙城斷面與藍(lán)天水閘接近，斷面水位過(guò)程與藍(lán)天水閘的閘前水位存在較大的相關(guān)性；帆游、塘下和海城斷面離下邊界水閘距離較遠(yuǎn)，需直接與各分區(qū)產(chǎn)匯流過(guò)程建立相關(guān)性求得。因此不同斷面的洪水位過(guò)程需采用不同的因子組合進(jìn)行模擬。

表1 溫瑞平原流域各斷面水位預(yù)報(bào)因子篩選結(jié)果

4.2 洪水位預(yù)報(bào)結(jié)果分析 基于篩選出的預(yù)報(bào)因子組合，構(gòu)建混合回歸模型。采用1999年9月3日—6日、2005年7月18日—21日、2007年9月18日—21日和2009年8月7日—10日共4場(chǎng)洪水作為率定期訓(xùn)練模型參數(shù)，對(duì)2010年7月24日—27日和2013年10月6日—9日對(duì)應(yīng)的9個(gè)斷面的洪水位進(jìn)行預(yù)報(bào)研究。通過(guò)對(duì)比發(fā)現(xiàn)，西山與水心殿斷面洪水過(guò)程相似。這是由于西山與水心殿斷面距離較近，同處于北排甌江的格局片區(qū)，河網(wǎng)水系連通，且洪水來(lái)源相同。永強(qiáng)和沙城斷面、帆游和塘下斷面也存在上述情況。由此可以得出，當(dāng)進(jìn)行沿海平原河網(wǎng)洪水預(yù)報(bào)工作時(shí)，可先結(jié)合洪水來(lái)源與排水走向，分析平原河網(wǎng)排澇格局，在此基礎(chǔ)上合理布置特征斷面，對(duì)其分別構(gòu)建混合回歸模型，即可快速模擬出整個(gè)平原的洪水位空間分布及洪水演進(jìn)情況。

為節(jié)省篇幅，本文分別從各片區(qū)中選擇一個(gè)代表站，以康華僑東、永強(qiáng)和塘下端面為例，對(duì)比六個(gè)場(chǎng)次洪水的實(shí)測(cè)與模擬過(guò)程，如圖3所示?？梢钥闯觯诼识ㄆ?，3個(gè)代表站水位的模擬值與實(shí)測(cè)值演變趨勢(shì)相同，且對(duì)水位峰值出現(xiàn)時(shí)間和峰值大小模擬效果較好，其中2009年模擬效果最好；但對(duì)低水位及退水段模擬不佳，尤其是1999年，低水位時(shí)模擬值偏大，水位發(fā)生波動(dòng)時(shí)的波谷值模擬也偏大，有待改進(jìn)。

圖3 溫瑞平原特征斷面洪水位實(shí)測(cè)與模擬對(duì)比過(guò)程圖

為了更進(jìn)一步評(píng)價(jià)模型的準(zhǔn)確性，本文統(tǒng)計(jì)了每個(gè)特征斷面每場(chǎng)洪水的特征值進(jìn)行精度評(píng)價(jià)，本研究選取的預(yù)報(bào)精度評(píng)價(jià)指標(biāo)為確定性系數(shù)（Dc）和平均合格率（QR），統(tǒng)計(jì)結(jié)果如表2所示?？梢钥闯?，洪峰水位計(jì)算值普遍比實(shí)際值偏大，但其差值小于0.10 m；峰現(xiàn)時(shí)間錯(cuò)時(shí)基本介于3 h以內(nèi)；確定性系數(shù)均大于0.70，部分洪水場(chǎng)次確定性系數(shù)甚至高達(dá)0.97；合格率基本高于70%，預(yù)報(bào)效果整體較好。

表2 混合回歸模型洪水位計(jì)算結(jié)果

針對(duì)溫瑞平原9個(gè)特征斷面率定期和檢驗(yàn)期分別進(jìn)行精度評(píng)價(jià)，評(píng)價(jià)結(jié)果詳見(jiàn)表3?？梢钥闯觯识ㄆ诘钠骄_定性系數(shù)為0.92，檢驗(yàn)期為0.82，均處于0.80以上；率定期的平均合格率達(dá)80.71%，檢驗(yàn)期平均合格率為75.60%，預(yù)報(bào)精度整體較高。從預(yù)報(bào)成果來(lái)說(shuō)，采用混合回歸模型進(jìn)行洪水位模擬滿足《水文情報(bào)預(yù)報(bào)規(guī)范》（GB/T 22482-2008）中規(guī)定的乙級(jí)預(yù)報(bào)精度標(biāo)準(zhǔn)，即合格率QR介于70.0%～85.0%，確定性系數(shù)Dc介于0.70～0.90，效果良好，可用于作業(yè)預(yù)報(bào)。

表3 混合回歸模型洪水位計(jì)算精度評(píng)價(jià)

值得注意的是，從圖3可以看出，康華僑東斷面洪水位過(guò)程并不像山區(qū)河道具有明顯的漲水段、退水段，而是波動(dòng)劇烈，除了明顯的峰值之外，還存在周期性的波動(dòng)情況。這是由于康華僑東斷面距離閘門很近，其水位受洪水與潮水雙重作用的影響。結(jié)合其他斷面水位過(guò)程得知，距離閘門越近，水位受閘門啟閉影響越大，波動(dòng)越劇烈。說(shuō)明在進(jìn)行沿海平原河網(wǎng)洪水預(yù)報(bào)研究時(shí)，潮位頂托和閘門控制是不可忽視的影響。

5 結(jié)論與建議

本文以溫瑞平原流域?yàn)槔?，采取物理成因與數(shù)理統(tǒng)計(jì)分析相結(jié)合的思路對(duì)沿海平原河網(wǎng)洪水位預(yù)報(bào)工作進(jìn)行研究，得到結(jié)論如下：

（1）根據(jù)平原河網(wǎng)洪水過(guò)程的時(shí)空變化特征，通過(guò)利用產(chǎn)匯流模型計(jì)算求得各水文分區(qū)洪水過(guò)程，與下游潮位水位過(guò)程一起形成預(yù)報(bào)因子庫(kù)，并在此基礎(chǔ)上針對(duì)各特征斷面構(gòu)建相應(yīng)的混合回歸模型。此方法思路清晰，計(jì)算簡(jiǎn)單，運(yùn)行穩(wěn)定，可為平原感潮河網(wǎng)洪水位預(yù)報(bào)提供一條有別于傳統(tǒng)洪水預(yù)報(bào)模式的快速建模途徑。

（2）利用2010年和2013年2場(chǎng)洪水的實(shí)測(cè)資料，對(duì)模型進(jìn)行驗(yàn)證。結(jié)果表明，驗(yàn)證洪水模擬值與實(shí)測(cè)值演變趨勢(shì)相同，且對(duì)水位峰值出現(xiàn)時(shí)間和峰值大小模擬效果較好，確定性系數(shù)為0.82，合格率為75.60%，預(yù)報(bào)精度較高。混合回歸模型在溫瑞平原河網(wǎng)洪水位預(yù)報(bào)中可以取得較好的效果。

（3）模型仍存在一定的優(yōu)化空間。例如，若河網(wǎng)中建有閘、堰等水工建筑物，在篩選因子和構(gòu)建模型時(shí)，還應(yīng)考慮水工建筑物對(duì)洪水的調(diào)蓄影響；在預(yù)報(bào)模型對(duì)低水位及退水段模擬不佳的問(wèn)題上，可嘗試將洪水過(guò)程分段考慮，分別針對(duì)高水位和低水位進(jìn)行預(yù)報(bào)因子篩選和混合回歸模型構(gòu)建工作。另外，混合回歸模型的預(yù)測(cè)精度與輸入的歷史數(shù)據(jù)量有關(guān)，歷史洪水樣本越多，模型參數(shù)越向最優(yōu)值靠近，因此在實(shí)際預(yù)報(bào)工作中，還需考慮實(shí)現(xiàn)模型的自我學(xué)習(xí)功能，這些有待日后的研究。