蒲寶卿

研究黃河水沙通量的變化規(guī)律對(duì)沿黃流域的環(huán)境治理、氣候變化和人民生活的影響，以及對(duì)優(yōu)化黃河流域水資源分配、協(xié)調(diào)人地關(guān)系、調(diào)水調(diào)沙、防洪減災(zāi)等方面都具有重要的理論指導(dǎo)意義.

目前許多水利及相關(guān)專家都高度關(guān)注黃河調(diào)水調(diào)沙問題，不斷研究探索如何降低黃河河床，改變黃河地上懸河的現(xiàn)象.邱學(xué)紹等［1］采用SAS軟件作線性回歸，確定了排沙量與水流量的函數(shù)關(guān)系，給出了SAS回歸模型.饒明貴等［2］依據(jù)灰色數(shù)學(xué)理論，從等時(shí)距GM（1，1）模型理論出發(fā)，建立了水流量和排沙量的非等時(shí)距GM（1，1）模型.張?jiān)刚碌龋?］運(yùn)用灰色系統(tǒng)理論，采用DGM、Verhu lst模型建立了流量與小浪底水庫排沙速度之間的關(guān)系，并運(yùn)用GM（1，1）殘差模型提高模擬精度.韓中庚［4］建立三次B樣條數(shù)學(xué)模型，可用于研究黃河水沙通量問題.曹明偉等［5］通過對(duì)黃河水流量和排沙量變化過程的分析得到函數(shù)關(guān)系，根據(jù)插值和擬合理論，利用Matlab軟件得到排沙量與時(shí)間、排沙量與流量的函數(shù)關(guān)系式，并進(jìn)行數(shù)值積分得到總排沙量.這些模型的建立與預(yù)測(cè)分析都為制定黃河調(diào)水調(diào)沙的實(shí)施政策提供了參考.

本文根據(jù)小浪底水庫下游黃河某水文站2016—2021年的水位、流量與含沙量等實(shí)際監(jiān)測(cè)數(shù)據(jù)［6］，設(shè)計(jì)實(shí)驗(yàn)、建立模型進(jìn)行預(yù)測(cè)，最終提出了水沙通量監(jiān)測(cè)的優(yōu)化策略.

1 實(shí)驗(yàn)設(shè)計(jì)

一是構(gòu)建含沙量隨時(shí)間變化圖，以及含沙量與水位、水流量關(guān)系的散點(diǎn)圖.在此基礎(chǔ)上確定含沙量與各因素之間的關(guān)系，通過相關(guān)性分析、回歸分析等進(jìn)行量化分析.水沙通量的突變性可以使用Mann-Kendall非參數(shù)檢驗(yàn)，確定水沙通量時(shí)間序列的突變年份.水沙通量時(shí)間序列的長(zhǎng)期趨勢(shì)、規(guī)律的季節(jié)項(xiàng)、不規(guī)律的周期項(xiàng)可以使用奇異譜分析（Singular Spectrum Analysis，SSA）進(jìn)行分解.

二是對(duì)季節(jié)項(xiàng)和周期項(xiàng)計(jì)算周期時(shí)間，將季節(jié)項(xiàng)、周期項(xiàng)與實(shí)際的氣候規(guī)律、政策因素等進(jìn)行關(guān)聯(lián)，以便更深入研究水沙通量的變化規(guī)律.對(duì)奇異譜分解的每個(gè)項(xiàng)分別進(jìn)行分析預(yù)測(cè)，在每個(gè)項(xiàng)規(guī)律明確的情況下，直接使用相應(yīng)的函數(shù)進(jìn)行擬合，在規(guī)律不明確的情況下，使用長(zhǎng)時(shí)間序列預(yù)測(cè)（Long-Term Series Forecasting，LTSF）進(jìn)行預(yù)測(cè).

三是根據(jù)調(diào)水調(diào)沙前后河底高程的變化規(guī)律分析小浪底水庫調(diào)水調(diào)沙的實(shí)際效果，進(jìn)而通過回歸分析構(gòu)建河底高程變化與排沙速度的關(guān)系.根據(jù)排沙速度的預(yù)測(cè)值和排沙速度與河底高程變化規(guī)律，預(yù)測(cè)未來10年河底高程的變化.

2 模型的建立與求解

2.1 含沙量與時(shí)間、水位、水流量的關(guān)系

①含沙量與時(shí)間的關(guān)系.因不同年份含沙量的平均值差異較大，為方便展示每年內(nèi)含沙量隨時(shí)間變化規(guī)律，首先繪制含沙量隨時(shí)間變化圖如圖1所示.從圖1可以看出，大部分年份含沙量數(shù)據(jù)在6—7月份間存在峰值，可能受小浪底水庫調(diào)水調(diào)沙政策的影響.2016—2021年年均含沙量分別為1.012 5 kg/m3、1.045 1 kg/m3、6.544 1 kg/m3、5.764 5 kg/m3、6.991 6 kg/m3、5.502 8 kg/m3.含沙量在2016—2017年相對(duì)較低，而在2018年顯著提高，并在2018—2021年期間保持相對(duì)平穩(wěn).

圖1 2016—2021年含沙量隨時(shí)間變化圖

②含沙量與水位的相關(guān)性分析及回歸分析.分別對(duì)2016—2021年各年份的含沙量與水位數(shù)據(jù)繪制散點(diǎn)圖，兩者的線性回歸分析和相關(guān)性分析分別如圖2及表1所示.

表1 含沙量與水位的相關(guān)系數(shù)

圖2 含沙量與水位散點(diǎn)圖及線性回歸分析

從圖2和表1可以看出，含沙量與水位存在較強(qiáng)的線性關(guān)系.2016—2020年，線性函數(shù)擬合優(yōu)度在0.40～0.62之間，但是2021年的擬合優(yōu)度僅為0.082 113.2016—2020年含沙量與水位相關(guān)系數(shù)均大于0.63，兩者存在顯著的正相關(guān)性.2021年含沙量與水位相關(guān)系數(shù)僅為0.292 744.

③含沙量與水流量的關(guān)系.分別對(duì)2016—2021年各年份的含沙量與水流量數(shù)據(jù)繪制散點(diǎn)圖，對(duì)兩者的線性回歸分析和相關(guān)性分析分別如圖3及表2所示.

表2 含沙量與水流量的相關(guān)系數(shù)

圖3 含沙量與水流量散點(diǎn)圖及線性回歸分析

從圖3和表2可以看出，水流量與含沙量的相關(guān)系數(shù)相對(duì)水位更高，其中2017年兩者的相關(guān)系數(shù)達(dá)到0.828 269.含沙量與水位的相關(guān)性可能是含沙量與水流量相關(guān)性的衍生［7］.

④總水流量和總排沙量估計(jì).含沙量的觀測(cè)數(shù)據(jù)相對(duì)水位和流量觀測(cè)數(shù)據(jù)存在大量缺失，為了更準(zhǔn)確地評(píng)估總排沙量，使用線性插值補(bǔ)全缺失的含沙量觀測(cè)數(shù)據(jù).另外，為了控制每天觀測(cè)次數(shù)的不同對(duì)總水流量和排沙量的影響，使用每天觀測(cè)值的均值代表當(dāng)天的觀測(cè)值.2016—2021年年總水流量和年總排沙量的估計(jì)值如表3所示.

表3 2016-2021年年總水流量和年總排沙量估計(jì)值

從表3可以看出，總水流量在2016—2021年期間總體呈現(xiàn)上升趨勢(shì)，2021年該值估計(jì)為4.741 8×1010m3.總排沙量在2020年達(dá)到峰值，數(shù)值為3.487 264×1011kg.

2.2 水沙通量的突變性、季節(jié)性和周期性

①水沙通量的Mann-Kendall突變性檢驗(yàn).通過Mann-Kendall非參數(shù)檢驗(yàn)法對(duì)年總水流量和年總排沙量的突變年份進(jìn)行確認(rèn)［8］，結(jié)果如圖4所示.從圖4可以看出，在年總水流量的Mann-Kendall統(tǒng)計(jì)量曲線中，在0.05置信度區(qū)間內(nèi)，UF和UB曲線在2017—2018年相交.年總排沙量在2017年出現(xiàn)交點(diǎn).因此，相應(yīng)的時(shí)間節(jié)點(diǎn)被推斷為水沙時(shí)間序列的突變年份.

圖4 年總水流量和年總排沙量的Mann-Kenda ll非參數(shù)檢驗(yàn)

②水沙通量的季節(jié)性及周期性的SSA.水沙通量的時(shí)間序列趨勢(shì)如圖5所示.為發(fā)現(xiàn)明顯的季節(jié)性因素，進(jìn)一步使用SSA對(duì)水流量和排沙量的時(shí)間序列進(jìn)行分解.SSA能夠根據(jù)觀測(cè)到的時(shí)間序列構(gòu)造軌跡矩陣，并對(duì)此進(jìn)行分解和重構(gòu)［9-10］，從而提取長(zhǎng)期趨勢(shì)因素、周期性因素、噪聲因素等不同成分.

圖5 水流量、排沙量趨勢(shì)的時(shí)間序列圖

使用Matlab軟件內(nèi)建函數(shù)trenddecomp對(duì)水流量和排沙量進(jìn)行SSA，水流量時(shí)間序列中發(fā)現(xiàn)3個(gè)季節(jié)性因素，其中方差最大的季節(jié)性因素的周期為1年，推斷與氣候條件相關(guān)，對(duì)應(yīng)于枯水季—豐水季循環(huán).排沙量時(shí)間序列中方差最大的季節(jié)性因素周期同樣為1年，與水流量時(shí)間序列一致.另外在殘差中發(fā)現(xiàn)1個(gè)周期性因素，2018—2021年每年的6—7月份間，水流量提升了約2 000 m3/s.這一周期性因素推測(cè)與小浪底水庫的調(diào)水調(diào)沙政策相關(guān).

2.3 水沙通量的變化趨勢(shì)及監(jiān)測(cè)方案

前文使用SSA對(duì)水流量和排沙量的時(shí)間序列進(jìn)行了分解，發(fā)現(xiàn)了水沙通量的長(zhǎng)期趨勢(shì)、與氣候相關(guān)的季節(jié)性因素、與小浪底水庫調(diào)水調(diào)沙相關(guān)的周期性因素［11］.下文將針對(duì)SSA中的長(zhǎng)期趨勢(shì)項(xiàng)、季節(jié)項(xiàng)、周期項(xiàng)分別進(jìn)行擬合預(yù)測(cè)，然后整合各項(xiàng)預(yù)測(cè)得到水沙通量的整體預(yù)測(cè)值［12-13］.由于各項(xiàng)的規(guī)律性非常明顯，結(jié)合長(zhǎng)短期記憶網(wǎng)絡(luò)（Long-Short Term Memory，LSTM）在較長(zhǎng)期時(shí)間下預(yù)測(cè)的不準(zhǔn)確性，本文使用相應(yīng)的函數(shù)直接進(jìn)行擬合.

①基于SSA的水流量預(yù)測(cè).對(duì)水流量SSA中解析的趨勢(shì)項(xiàng)、季節(jié)項(xiàng)、周期項(xiàng)等分別進(jìn)行擬合.其中長(zhǎng)期趨勢(shì)項(xiàng)用線性函數(shù)擬合，擬合函數(shù)為：

式中：x表示時(shí)間，y表示水流量，k和b為待擬合參數(shù).

季節(jié)項(xiàng)使用正弦函數(shù)與非平穩(wěn)項(xiàng)聯(lián)合預(yù)測(cè)，擬合函數(shù)為：

式中：k1，k2，k3，ω為待擬合參數(shù).

調(diào)水調(diào)沙的周期項(xiàng)通過以下函數(shù)定義：

水流量時(shí)間序列中長(zhǎng)期趨勢(shì)項(xiàng)和季節(jié)項(xiàng)的擬合結(jié)果如圖6所示.

圖6 水流量時(shí)間序列中長(zhǎng)期趨勢(shì)項(xiàng)和季節(jié)項(xiàng)擬合圖

通過擬合最終得到的預(yù)測(cè)函數(shù)為：

式中：F（x）表示調(diào)水調(diào)沙相關(guān)的周期項(xiàng)，x表示距離2016年1月1日的天數(shù)，y的單位為m3/s.

基于SSA的未來兩年水流量預(yù)測(cè)如圖7所示.

圖7 基于SSA的未來兩年水流量預(yù)測(cè)圖

從圖7可以看出，水流量在每年的7月中旬和10月上旬存在兩個(gè)顯著的高峰期.2022年的高峰期預(yù)測(cè)發(fā)生在2022年7月12日和2022年10月2日.2023年的高峰期預(yù)測(cè)發(fā)生在2023年7月8日和2023年10月2日.

②基于SSA的排沙量預(yù)測(cè).與排水量預(yù)測(cè)過程相似，分別對(duì)SSA中各成分進(jìn)行預(yù)測(cè)并整合，最終得到排沙量的預(yù)測(cè)函數(shù)為：

式中：F（x）表示調(diào)水調(diào)沙相關(guān)的周期項(xiàng)，x表示距離2016年1月1日的天數(shù)，y的單位為kg/s.

③采樣監(jiān)測(cè)方案設(shè)計(jì).根據(jù)預(yù)測(cè)結(jié)果，未來兩年排水量的峰值出現(xiàn)在7月中旬，排沙量的峰值出現(xiàn)在7月中下旬，另外排水量在每年4月份和1月份存在兩個(gè)小高峰.因此，為及時(shí)掌握水沙通量動(dòng)態(tài)變化，又能最大程度減少監(jiān)測(cè)成本，應(yīng)在接近預(yù)估峰值的時(shí)間段提高監(jiān)測(cè)頻率［14］.另外，可以根據(jù)監(jiān)測(cè)數(shù)據(jù)的累積對(duì)水沙通量進(jìn)行實(shí)時(shí)預(yù)測(cè)，動(dòng)態(tài)設(shè)定監(jiān)測(cè)頻率.

2.4 小浪底水庫調(diào)水調(diào)沙實(shí)際效果分析

由于2016—2017年6—7月份間水流量增加不顯著，對(duì)2019年及以后的相關(guān)數(shù)據(jù)進(jìn)行分析.觀測(cè)數(shù)據(jù)中給出了不同采樣時(shí)間下起點(diǎn)距離與河底高程的關(guān)系，不同時(shí)間下河底高程比較如圖8所示.

圖8 河底高程比較圖

從圖8可以看出，2019年4月13日、2019年10月15日、2020年3月19日及2021年3月14日的平均河底高程分別為44.969 1 m、43.864 4 m、45.089 8 m、45.023 4 m.因?yàn)椴糠制瘘c(diǎn)距離的數(shù)據(jù)缺失，僅計(jì)算數(shù)據(jù)存在區(qū)段的平均河底高程并繪圖.平均河底高程在2019年4月13日至2019年10月15日期間下降了1.104 7 m，期間經(jīng)歷了6—7月份小浪底水庫調(diào)水調(diào)沙期，說明調(diào)水調(diào)沙對(duì)清理河底淤沙有明顯效果.而2019年10月15日至2020年3月19日期間無調(diào)水調(diào)沙，平均河底高程由43.864 4 m上升到45.089 8 m，上升了1.225 4米，說明若不經(jīng)調(diào)水調(diào)沙，河底淤沙會(huì)快速上升.在2020年3月19日至2021年3月14日大約一年的時(shí)間內(nèi)，平均河底高程幾乎無變化，說明6—7月份調(diào)水調(diào)沙能夠保障河底高程處于動(dòng)態(tài)平衡中.

為了研究不進(jìn)行調(diào)水調(diào)沙下10年后河底高程預(yù)測(cè)，直接建立排沙速度與河底高程速度變化的聯(lián)系.根據(jù)表3中的數(shù)據(jù)和水沙通量的計(jì)算結(jié)果，繪制河底高程增長(zhǎng)速度與排沙速度的關(guān)系，如圖9所示.兩者之間有非常顯著的線性相關(guān)性.

圖9 排沙速度與河底高程增長(zhǎng)速度圖

經(jīng)線性方程擬合得到，排沙速度與河底高程增長(zhǎng)速度的關(guān)系為：

式中：x表示排沙速度（t/s），y表示河底高程增長(zhǎng)速度（m/d）.

當(dāng)排沙速度為10.519 5 t/s時(shí)，河底平均高程將保持不變.排沙速度小于10.519 5 t/s時(shí)，河底高程將增加；反之，河底高程將減小.

為預(yù)測(cè)不進(jìn)行調(diào)水調(diào)沙下，10年后河底高程變化，去掉6—7月份的調(diào)水調(diào)沙代表的周期項(xiàng)后，對(duì)排沙速度進(jìn)行預(yù)測(cè)，得到未來10年平均排沙速度為7.964 9 t/s時(shí)，河底高程平均增長(zhǎng)速度為0.002 2 m/d，10年后河底平均高程將增加8.03 m.

3 結(jié)語

根據(jù)黃河某水文站6年的相關(guān)監(jiān)測(cè)數(shù)據(jù)，建立含沙量與時(shí)間、水位、水流量的關(guān)系，估算這6年的年總水流量和年總排沙量.分析水沙通量的突變性、季節(jié)性、周期性等變化規(guī)律，構(gòu)建模型對(duì)水沙通量進(jìn)行量化分析，使用SSA法對(duì)水沙通量時(shí)間序列進(jìn)行了分解，解析出其中的長(zhǎng)期趨勢(shì)、與氣候相關(guān)的季節(jié)性規(guī)律和與調(diào)水調(diào)沙相關(guān)的周期性規(guī)律.

根據(jù)結(jié)果對(duì)未來兩年水沙通量進(jìn)行預(yù)測(cè)，在此基礎(chǔ)上提出水沙通量監(jiān)測(cè)的優(yōu)化策略.最后量化河底高程隨水沙通量的變化規(guī)律，分析了小浪底水庫調(diào)水調(diào)沙的實(shí)際效果，可為黃河水沙監(jiān)測(cè)計(jì)劃及調(diào)水調(diào)沙政策的制定提供參考.