黃河水沙監(jiān)測數(shù)據(jù)分析與模型構(gòu)建

蒲寶卿

研究黃河水沙通量的變化規(guī)律對沿黃流域的環(huán)境治理、氣候變化和人民生活的影響，以及對優(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é)理論，從等時距GM（1，1）模型理論出發(fā)，建立了水流量和排沙量的非等時距GM（1，1）模型.張愿章等［3］運用灰色系統(tǒng)理論，采用DGM、Verhu lst模型建立了流量與小浪底水庫排沙速度之間的關(guān)系，并運用GM（1，1）殘差模型提高模擬精度.韓中庚［4］建立三次B樣條數(shù)學(xué)模型，可用于研究黃河水沙通量問題.曹明偉等［5］通過對黃河水流量和排沙量變化過程的分析得到函數(shù)關(guān)系，根據(jù)插值和擬合理論，利用Matlab軟件得到排沙量與時間、排沙量與流量的函數(shù)關(guān)系式，并進(jìn)行數(shù)值積分得到總排沙量.這些模型的建立與預(yù)測分析都為制定黃河調(diào)水調(diào)沙的實施政策提供了參考.

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

1 實驗設(shè)計

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

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

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

2 模型的建立與求解

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

①含沙量與時間的關(guān)系.因不同年份含沙量的平均值差異較大，為方便展示每年內(nèi)含沙量隨時間變化規(guī)律，首先繪制含沙量隨時間變化圖如圖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年相對較低，而在2018年顯著提高，并在2018—2021年期間保持相對平穩(wěn).

圖1 2016—2021年含沙量隨時間變化圖

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

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

圖2 含沙量與水位散點圖及線性回歸分析

從圖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)系.分別對2016—2021年各年份的含沙量與水流量數(shù)據(jù)繪制散點圖，對兩者的線性回歸分析和相關(guān)性分析分別如圖3及表2所示.

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

圖3 含沙量與水流量散點圖及線性回歸分析

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

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

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

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

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

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

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

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

圖5 水流量、排沙量趨勢的時間序列圖

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

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

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

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

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

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

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

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

水流量時間序列中長期趨勢項和季節(jié)項的擬合結(jié)果如圖6所示.

圖6 水流量時間序列中長期趨勢項和季節(jié)項擬合圖

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

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

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

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

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

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

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

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

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

由于2016—2017年6—7月份間水流量增加不顯著，對2019年及以后的相關(guān)數(shù)據(jù)進(jìn)行分析.觀測數(shù)據(jù)中給出了不同采樣時間下起點距離與河底高程的關(guān)系，不同時間下河底高程比較如圖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.因為部分起點距離的數(shù)據(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)沙對清理河底淤沙有明顯效果.而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)沙，河底淤沙會快速上升.在2020年3月19日至2021年3月14日大約一年的時間內(nèi)，平均河底高程幾乎無變化，說明6—7月份調(diào)水調(diào)沙能夠保障河底高程處于動態(tài)平衡中.

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

圖9 排沙速度與河底高程增長速度圖

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

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

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

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

3 結(jié)語

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

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