劉濤 周曼 胡挺 紀國良
摘要:三峽水庫蓄水以來,由于入庫流量統(tǒng)計和地形資料測量等存在系統(tǒng)性誤差,通過求解圣維南方程組計算的長壽站水位精度有時難以滿足實際應用需求。對此,采用分段套索最小角回歸交叉驗證(Lasso Least angle regression Cross Validation,LassoLarsCV)算法預測長壽站水位。首先收集了2009~2019年相關水位和流量數(shù)據(jù)并進行標準化處理,然后根據(jù)壩前水位和入庫流量對數(shù)據(jù)進行分類;最后針對每一個類別,從2009~2018年的數(shù)據(jù)中隨機抽取80%用于訓練線性模型,剩余20%數(shù)據(jù)用于檢驗并調(diào)整分段方案,以2019年的數(shù)據(jù)進行測試,采用LassoLarsCV算法模型求解。結(jié)果表明:采用LassoLarsl算法得到的長壽站水位誤差可以控制在1.2 m以內(nèi),99.47%的水位誤差在0.5 m以內(nèi),效果優(yōu)于傳統(tǒng)水動力學方法。LassoLarsCV可直接研究出、入庫流量和壩前水位與長壽站水位的映射關系,降低了模型對輸入數(shù)據(jù)準確性的要求,同時提高了長壽站水位的計算精度。
關 鍵 詞:
水位預測; LassoLarsCV算法; 圣維南方程組; 線性模型; 水動力學; 三峽庫區(qū)
中圖法分類號: P332.3
文獻標志碼: A
DOI:10.16232/j.cnki.1001-4179.2021.05.010
1 研究背景
三峽水庫蓄水以來,根據(jù)長期積累的研究成果和實踐經(jīng)驗,發(fā)現(xiàn)長壽段是回水淹沒的敏感區(qū)域。該區(qū)域長壽站水位是否超過土地征用線或移民遷移線,往往代表著庫尾是否發(fā)生淹沒,因此需要重點關注。在以往的研究中,多采用洪水演進模型計算長壽站水位。但是由于三峽水庫覆蓋地域廣、河網(wǎng)復雜,模型的輸入邊界難以準確獲取,長壽站水位的計算精度有時難以滿足實際應用需求。因此,有必要進一步研究高精度的長壽站水位預測方法。
目前,已有許多研究致力于提高洪水演進模型的計算精度。長江科學院在計算三峽水庫洪水傳播過程中,一方面增加了水庫斷面測量點數(shù)量,降低了斷面概化誤差;另一方面不斷提升入庫流量的預報水平,并將更多支流加入到模型中。在糙率率定方面,楊世孝等[1]提出了一種反求糙率的方式;陳素紅等[2]建立了基于多親遺傳算法的河道糙率率定模型;陳一帆等[3]以糙率和水力狀態(tài)量作為河網(wǎng)非線性動態(tài)系統(tǒng)狀態(tài)變量,采用擴展卡爾曼濾波構(gòu)建了結(jié)合糙率動態(tài)校正的河網(wǎng)水情數(shù)據(jù)同化模型。這些研究已經(jīng)取得了較大的進展,但它們?nèi)匀粚吔巛斎霔l件的準確性有較高要求,對于三峽水庫這類大型河網(wǎng),基本上都存在2 m左右的誤差,其在普適性和計算精度方面仍有待進一步提高。
由Robert Tibshirani提出的Lasso方法發(fā)展而來的套索最小角回歸交叉驗證(Lasso Least angle regression Cross Validation,LassoLarsCV)算法兼顧了計算速度的高效性和擬合結(jié)果的精確性,是目前應用最廣泛的一種線性學習算法。因此本文嘗試采用分段LassoLarsCV算法計算長壽站水位,該算法直接學習出、入庫流量和壩前水位到長壽站水位的映射關系,只需要當前出、入庫流量的統(tǒng)計誤差與歷史數(shù)據(jù)的誤差保持在一定范圍內(nèi)即可,對數(shù)據(jù)準確性的要求較洪水演進模型低。
LassoLarsCV算法采用最小角回歸(Least angle regression)法求解式(2) 并進行交叉驗證。最小角回歸法是一種快速進行特征選擇和回歸系數(shù)計算的迭代算法,其核心思想是將回歸目標向量依次分解為若干組特征向量的線性組合,最終使得與所有特征均線性無關的殘差向量最小[8]。具體來說,算法首先找到與因變量y相關度最高的xk(余弦距離最大),在已選擇變量的solution path上前進直到有另一個變量xt,使得這兩個變量的所有取值實例構(gòu)成的特征向量Xk、Xt與當前殘差的相關度是一樣的。 然后重復這個過程,直到殘差向量足夠得小,或者說所有的變量都已經(jīng)取完,算法停止。目前,最小角回歸法在數(shù)據(jù)分析中十分流行,比如用于光譜分析模型的建立、結(jié)構(gòu)地震需求重要性度量分析、全基因組選擇等研究工作中[9-11]。
LassoLarsCV算法的建模和求解可以通過調(diào)用Python中sklearn包linear_model模塊中的LassoLarsCV類進行[12]。
3 長壽站水位預測
影響長壽站水位的水文因素眾多,以這些因素作為自變量,則長壽站水位是它們的復雜函數(shù)。理論上可通過對相關變量取值范圍的逐漸細分,在每個局部采用線性回歸模型擬合該復雜函數(shù)。在一個局部細分范圍內(nèi),樣本數(shù)量越充足、分布越均勻,LassoLarsCV算法求解得到的線性回歸模型擬合預測的效果越好。
3.1 變量選擇
根據(jù)水動力學圣維南方程組原理,在已知水庫干、支流的入庫流量和壩前水位變化過程條件下,可以計算出庫區(qū)內(nèi)所有斷面的流量和水位過程[13]。在計算過程中,不同斷面的流量和水位之間存在相互關聯(lián)作用不是各自獨立的。因此,在使用LassoLarsCV算法來研究長壽站水位的預測問題時,需要找到關鍵的影響因素。
三峽庫區(qū)主要水文站點從上游到下游的順序依次為朱沱站、北碚站、寸灘站、長壽站、武隆站、三峽水庫中心站,其中朱沱站、寸灘站、長壽站、三峽水庫中心站為干流站點,北碚站、武隆站為支流站點,具體如圖1所示。因此從長壽站所處的位置來看,寸灘站水文信息可以比較全面反映來自長壽上游的朱沱站和北碚站(兩個主要站點)的水文因素對長壽站水位的綜合影響;同時,長壽站下游的武隆站和三峽水庫中心站的相關水文信息反映了下游對長壽站水位的頂托效應。將上述理論分析結(jié)合實踐,初步確定影響長壽站水位的重要因素為寸灘流量、武隆流量、三峽水庫出庫流量、三峽水庫平滑流量(一種可以通過預報得到的入庫流量)、鳳凰山水位(由三峽水庫中心站給出的壩前水位,位于鳳凰山而得名)。
在具體的模型變量選擇時,不僅要考慮提高長壽站水位預測精度,而且要考慮實際預報的能力。當前三峽水庫入庫流量預報的時間間隔為6 h,兩個時間間隔之間的入庫流量和水位均未知,插值法可能在一般情況下是適用的,但是在汛期大流量情況下,1 h內(nèi)流量和水位變化劇烈,存在很大的不確定性。因此上述各因素的取值只能取6 h間隔。
在確定影響長壽站水位的主要影響因素以及取值時間間隔后,還需要選擇每個因素數(shù)據(jù)的時段數(shù)。這里“時段數(shù)”是指在預測長壽站水位時,每個因素需要使用的時間點個數(shù),例如當時段數(shù)為2時,假設要對t時刻長壽站水位進行預測,各主要影響因素的取值為t-6時刻和t時刻的值??紤]到既能反映水文過程又不造成變量冗余,選擇兩個時段是比較合適的。
綜上所述,長壽站水位預測的分段LassoLarsCV算法的特征變量為t-6時刻以及t時刻的寸灘流量、武隆流量、三峽水庫出庫流量、三峽水庫平滑流量、鳳凰山水位,共計10個變量,標簽為t時刻的實測長壽站水位。
3.2 數(shù)據(jù)處理
在長壽站水位預測問題中,存在流量和水位兩種數(shù)值,它們之間量級相差較大。為了避免流量因子被LassoLarsCV算法錯誤地壓縮掉,在將樣本帶入模型中進行計算時會先進行數(shù)據(jù)標準化的預處理,在這里采用的標準化處理方法具體為
標準化水位數(shù)據(jù)=原水位數(shù)據(jù)-原最低水位原最高水位-原最低水位
標準化流量數(shù)據(jù)=原流量數(shù)據(jù)-原最小流量原最大流量-原最小流量
根據(jù)三峽水庫的運行和數(shù)據(jù)積累情況,選取2009年1月1日00:00至2019年11月18日23:00的數(shù)據(jù)進行擬合訓練和預測。其中2009年1月1日00:00至2018年12月31日23:00共87 648條數(shù)據(jù),形成的所有樣本按照t時刻鳳凰山水位和t時刻三峽水庫平滑流量兩個特征變量取值范圍進行細分,細分的具體操作是將當前時刻鳳凰山水位和三峽水庫平滑流量兩個特征變量分別從樣本中各自最小值開始,以固定值Δz,Δq為間隔進行網(wǎng)格劃分,滾動完成對所有樣本的分類。然后在每個分類下隨機抽取80%的樣本數(shù)據(jù)用于學習訓練,剩下20%用于驗證,根據(jù)其預測效果反饋調(diào)整Δz,Δq。2019年1月1日00:00至2019年11月18日23:00共7 728條數(shù)據(jù),全部用于測試模型的預測效果。
3.3 模型評價指標
為了合理評價通過訓練得到的線性回歸模型的性能,考慮通過以下3個指標進行量化展示。
3.4 結(jié)果分析
經(jīng)過計算,找到了一個比較合適的分類間隔Δz=1 m,Δq=20 000 m3/s,即t時刻鳳凰山水位每隔1 m且三峽水庫平滑流量每隔20 000 m3/s一個分類。
3.4.1 擬合和預測效果
表1是分段LassoLarsCV算法得到的線性回歸模型分別在訓練集、驗證集和測試集上對長壽站水位進行預測的誤差區(qū)間統(tǒng)計情況,可以看出最大誤差在1.20 m以內(nèi)。經(jīng)過計算,模型的平均絕對誤差為0.10 m,均方誤差為0.017 5,各個分段的決定系數(shù)均在0.88以上。2019年的測試結(jié)果顯示:預測誤差在[-0.81 m,0.74 m]區(qū)間范圍內(nèi),在該年度最大來流(40 000~50 000 m3/s)期間預測誤差控制在0.45 m以內(nèi),效果較好。
圖2是2009~2019年長壽站水位擬合預測誤差分布,表2是2009~2019年長壽站水位擬合預測誤差絕對值分段統(tǒng)計。從11 a的情況來看,誤差主要集中在±0.5 m以內(nèi)。在所有的擬合預測結(jié)果中,誤差絕對值小于0.5 m的數(shù)據(jù)比例高達99.47%,誤差絕對值大于或者等于0.5 m的數(shù)據(jù)僅占比0.53%。
3.4.2 LassoLarsCV算法與傳統(tǒng)水動力學法比較
在前面提到過傳統(tǒng)水動力學法在計算長壽站水位時最大誤差高達2 m,為了更加充分地說明這一點,在這一節(jié)將具體展示部分傳統(tǒng)水動力學法計算結(jié)果與本文結(jié)果的比較。2012,2014年和2017年是長江流域典型洪水年,均出現(xiàn)過較大流量或者較高水位,因此有必要用這3 a數(shù)據(jù)將傳統(tǒng)水動力學法和分段LassoLarsCV算法的結(jié)果與實測水位的結(jié)果進行對比,具體情況如圖3~5中所示。從圖3~5中不難看出,與本文采用的分段LassoLarsCV算法相比,傳統(tǒng)水動力學法的誤差分布范圍較大一些,而且振蕩較為劇烈、太不穩(wěn)定,最大誤差超過了2 m,效果存在一定差距。
3.4.3 零界點處誤差分析
本節(jié)重點分析長壽站實際水位高于175.00 m(在土地征用線附近)時的預測效果。對歷史數(shù)據(jù)進行統(tǒng)計分析,發(fā)現(xiàn)長壽站實際水位高于175.00 m主要有兩種情況:① 在入庫大流量(q>40 000 m3/s)時發(fā)生,主要集中在汛期附近;② 在鳳凰山高水位(z>173.00 m)時發(fā)生,主要集中在蓄水期末期附近。下面分別就這兩種情況進行說明。
2012年7月24日和2014年9月19日均是在入庫流量較大時長壽站水位高于175.70 m(土地征用線高程),經(jīng)過統(tǒng)計這兩個時段的預測誤差,可以得出以下結(jié)論:
在長壽站水位高于175.00 m時兩時段最大誤差分別為0.21 m和-0.07 m,最大誤差發(fā)生的時間均在水位從高位下降的過程中,整體誤差范圍分別為[-0.08 m,0.21 m]和[-0.07 m,0.06 m]。因此,長壽站水位高于175.00 m時,水位預測整體誤差都較小,當然這跟歷史樣本較少也有很大的關系,不足以說明用來預測超出歷史范圍的情況時也能有較好的效果。
預測結(jié)果在大多數(shù)時候存在滯后效應,即在長壽站實際水位上漲時的大多數(shù)情況下,預測的長壽站水位上漲相對較慢,誤差呈現(xiàn)為負值;而在長壽站實際水位下降時的大多數(shù)情況下,預測的長壽站水位下降相對較慢,誤差呈現(xiàn)為正值。這說明在長壽站水位上漲的過程中可以在預測值的基礎上加上一個正值進行修正。
當鳳凰山水位較高時,三峽水庫平滑流量通常較低,在此種情況下,當長壽站水位超過175.70 m時,2011,2014年和2017年的預測誤差最大值分為別-0.12,-0.37,-0.23 m,均在可接受范圍內(nèi)。
3.4.4 LassoLarsCV算法與深度學習方法、傳統(tǒng)水動力學法的區(qū)別
與深度學習方法相比,分段LassoLarsCV算法原理更易理解、模型更直觀、硬件要求低、訓練速度快。一方面,因為深度學習方法需整體進行訓練擬合,不同局部之間彼此聯(lián)系、相互影響,而本文的方法由于需要分段使得不同局部之間彼此獨立、缺乏聯(lián)系,不能相互提供有效的信息,所以在樣本稀少的區(qū)域,前者泛化預測效果較好,后者過擬合情況較為嚴重、泛化預測效果較差,而在超出歷史范圍較大的情況下兩者均不具備有效的預測能力。另一方面,正因為深度學習方法是從全局出發(fā)考慮問題,會在局部做出讓步以保持全局的擬合預測效果,所以在樣本數(shù)量較為充足、分布較為均勻的局部情況下本文方法的擬合預測效果要比深度學習方法稍好一些。
與傳統(tǒng)水動力學方法相比,分段LassoLarsCV算法充分挖掘利用了歷史數(shù)據(jù)提供的信息,需要的數(shù)據(jù)源少,計算過程簡便快捷。一方面,因為分段LassoLarsCV算法是從數(shù)據(jù)分析的角度處理問題,降低了對輸入數(shù)據(jù)準確度的要求,只需要預測和訓練時使用的數(shù)據(jù)屬于同一來源即可,在歷史樣本數(shù)量較為充足的情況下預測精度明顯提高。另一方面,在超出歷史范圍較大的情況下,不依賴于歷史數(shù)據(jù)提供信息的傳統(tǒng)水動力學方法又有其獨特的優(yōu)勢,這一點是分段LassoLarsCV算法和深度學習方法都無法媲美的。
4 結(jié) 論
針對長壽站水位預測問題,本文使用分段LassoLarsCV算法在每一個局部用線性模型去擬合逼近真實情況。計算結(jié)果顯示:最大誤差可以控制在±1.2 m以內(nèi),精度達到96%,且99%的情況下可以將誤差控制在±0.5 m以內(nèi);特別在高水位不小于170.00 m的情況下,即在長壽站水位最有可能超過土地征用線、移民遷移線時,誤差可以控制在[-0.37 m,0.33 m]區(qū)間內(nèi)。通過與深度學習方法和傳統(tǒng)水動力學方法的分析對比發(fā)現(xiàn):分段LassoLarsCV算法在歷史樣本數(shù)量較為充足、分布較為均勻的局部情況下擬合預測精度較高,因此具有實時更新知識庫的獨特效用;分段LassoLarsCV算法在歷史樣本數(shù)量相對較少的局部情況下,泛化預測效果不及模型更復雜、計算量更大的深度學習方法;在超出歷史范圍較大的水情下,包括分段LassoLarsCV算法和深度學習方法在內(nèi)的基于數(shù)據(jù)分析的預測技術仍無法取代傳統(tǒng)水動力學方法,因此在實踐應用中需要根據(jù)具體情況結(jié)合使用3種方法。
參考文獻:
[1] 楊世孝,肖子良.反求糙率的一種數(shù)值方法[J].數(shù)值計算與計算機應用,1994(4):247-260.
[2] 陳素紅,劉孟凱,邢領航.基于多親遺傳算法的渠道非恒定流糙率率定模型及其應用[J].水電能源科學,2013,31(5):81-83.
[3] 陳一帆,程海洋,萬曉麗,等.結(jié)合糙率校正的河網(wǎng)水情數(shù)據(jù)同化[J].水科學進展,2015,26(5):731-738.
[4] HOERL A E,KENNARD R W.Ridge Regression:Applications to Nonorthogonal Problems[J].Technometrics,1970,12(1);69-72.
[5] 王俊迪,許蘊山,彭芳,等.基于嶺回歸的紅外協(xié)同定位優(yōu)化算法[J].北京航空航天大學學報,2020,46(3):563-570.
[6] ROBERT T.Regression shrinkage and selection via the lasso [J].Journal of the Royal Statistical Society:Series B (Statistical Methodology),1996,58(1):267-288.
[7] 李小虹,常睿春.基于Lasso回歸模型的區(qū)域可持續(xù)發(fā)展研究[J].科技經(jīng)濟導刊,2020(26):234-235.
[8] RYAN T,ROB T,JONATHAN T,et al.Least angle regression[J].Annals of Statistics,2004,32(2):407-451.
[9] 熊芩,張若秋,李輝,等.最小角回歸算法(LAR)結(jié)合采樣誤差分布分析(SEPA)建立穩(wěn)健的近紅外光譜分析模型[J].分析測試學報,2018,37(7):778-783.
[10] 王秀振,錢永久,宋帥.基于隨機森林和最小角回歸的結(jié)構(gòu)地震需求重要性度量分析[J].振動與沖擊,2019,38(4):115-120.
[11] 孫嘉利,吳清太,溫陽俊,等.基于FASTmrEMMA、最小角回歸和隨機森林的全基因組選擇方法探索[J/OL].南京農(nóng)業(yè)大學學報(2020-08-18) [2020-10-19].http://kns.cnki.net/kcms/detail/32.1148.S.20200818.1418.002.html.
[12] SANNER M F.Python:a programming language for software integration and development[J].Journal of Molecular Graphics & Modelling,1999,17(1):57-61.
[13] 雒文生.河流水文學[M].北京:水利電力出版社,1992.
(編輯:江 文)