黃河水沙監測數據分析與模型構建

蒲寶卿

研究黃河水沙通量的變化規律對沿黃流域的環境治理、氣候變化和人民生活的影響，以及對優化黃河流域水資源分配、協調人地關系、調水調沙、防洪減災等方面都具有重要的理論指導意義.

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

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

1 實驗設計

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

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

三是根據調水調沙前后河底高程的變化規律分析小浪底水庫調水調沙的實際效果，進而通過回歸分析構建河底高程變化與排沙速度的關系.根據排沙速度的預測值和排沙速度與河底高程變化規律，預測未來10年河底高程的變化.

2 模型的建立與求解

2.1 含沙量與時間、水位、水流量的關系

①含沙量與時間的關系.因不同年份含沙量的平均值差異較大，為方便展示每年內含沙量隨時間變化規律，首先繪制含沙量隨時間變化圖如圖1所示.從圖1可以看出，大部分年份含沙量數據在6—7月份間存在峰值，可能受小浪底水庫調水調沙政策的影響.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年期間保持相對平穩.

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

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

表1 含沙量與水位的相關系數

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

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

③含沙量與水流量的關系.分別對2016—2021年各年份的含沙量與水流量數據繪制散點圖，對兩者的線性回歸分析和相關性分析分別如圖3及表2所示.

表2 含沙量與水流量的相關系數

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

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

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

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

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

2.2 水沙通量的突變性、季節性和周期性

①水沙通量的Mann-Kendall突變性檢驗.通過Mann-Kendall非參數檢驗法對年總水流量和年總排沙量的突變年份進行確認［8］，結果如圖4所示.從圖4可以看出，在年總水流量的Mann-Kendall統計量曲線中，在0.05置信度區間內，UF和UB曲線在2017—2018年相交.年總排沙量在2017年出現交點.因此，相應的時間節點被推斷為水沙時間序列的突變年份.

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

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

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

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

2.3 水沙通量的變化趨勢及監測方案

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

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

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

季節項使用正弦函數與非平穩項聯合預測，擬合函數為：

式中：k1，k2，k3，ω為待擬合參數.

調水調沙的周期項通過以下函數定義：

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

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

通過擬合最終得到的預測函數為：

式中：F（x）表示調水調沙相關的周期項，x表示距離2016年1月1日的天數，y的單位為m3/s.

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

圖7 基于SSA的未來兩年水流量預測圖

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

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

式中：F（x）表示調水調沙相關的周期項，x表示距離2016年1月1日的天數，y的單位為kg/s.

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

2.4 小浪底水庫調水調沙實際效果分析

由于2016—2017年6—7月份間水流量增加不顯著，對2019年及以后的相關數據進行分析.觀測數據中給出了不同采樣時間下起點距離與河底高程的關系，不同時間下河底高程比較如圖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.因為部分起點距離的數據缺失，僅計算數據存在區段的平均河底高程并繪圖.平均河底高程在2019年4月13日至2019年10月15日期間下降了1.104 7 m，期間經歷了6—7月份小浪底水庫調水調沙期，說明調水調沙對清理河底淤沙有明顯效果.而2019年10月15日至2020年3月19日期間無調水調沙，平均河底高程由43.864 4 m上升到45.089 8 m，上升了1.225 4米，說明若不經調水調沙，河底淤沙會快速上升.在2020年3月19日至2021年3月14日大約一年的時間內，平均河底高程幾乎無變化，說明6—7月份調水調沙能夠保障河底高程處于動態平衡中.

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

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

經線性方程擬合得到，排沙速度與河底高程增長速度的關系為：

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

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

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

3 結語

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

根據結果對未來兩年水沙通量進行預測，在此基礎上提出水沙通量監測的優化策略.最后量化河底高程隨水沙通量的變化規律，分析了小浪底水庫調水調沙的實際效果，可為黃河水沙監測計劃及調水調沙政策的制定提供參考.