大沽河水沙演化規律及驅動因素研究

胡萌　宋委　劉克金　吳振　雷炳霄　王毅

［關鍵詞］Mann-Kendall檢驗；突變檢驗；小波分析；氣候變化；徑流；泥沙；大沽河

［摘要］基于1956—2016年長系列水沙資料，利用MK、滑動T檢驗、累積距平曲線、小波分析等方法分析了大沽河水沙變化規律，在此基礎上，計算了氣候變化和人類活動對徑流量、輸沙量變化影響的貢獻率。結果表明：大沽河徑流量、輸沙量呈顯著減少趨勢，降水量呈不顯著減少趨勢；徑流量在1976年、輸沙量在1965年和1979年均發生了由多到少的突變；天然徑流量存在8 a、12 a、21 a三個變化周期，其中21 a是第一主周期；輸沙量周期變化存在5～8 a、9～17 a和18～30 a三類較明顯的時間尺度，變化的第一主周期是12 a；人類活動是造成徑流量、輸沙量減少的主要原因，人類活動對徑流量變化的貢獻率為76.8%～97.4%，對輸沙量變化的貢獻率為76.2%～98.4%。

［中圖分類號］ TV213.9［文獻標識碼］ A［文章編號］ 1000－0941（2023）08－0032－06

流域系統是對氣候變化與人類活動異常敏感的動態系統，水沙演化是其中最為活躍的部分[1]。徑流量和輸沙量的變化不僅關系到河流自身的發展演變，也反映了流域的環境狀況、水土流失程度及人類活動的影響[2-3]。

已有的研究主要關注大型流域系統水沙運移規律及其影響因素，而中小型流域系統的重要性往往被忽視或低估［4-6］。大沽河是膠東半島最大的河流，地理位置與氣候環境獨特，流域面積占青島市域總面積的42%，被譽為青島市的“母親河”。

目前，對大沽河水沙研究多集中于單一要素方面［7-10］，但尚未見針對大沽河流域開展氣候變化與人類活動雙重驅動下水沙耦合演化方面的研究。本研究以大沽河流域降水、徑流、泥沙為研究對象，分析水沙演化規律，定量研究氣候變化和人類活動對徑流量和輸沙量變化影響的貢獻率，以期為流域環境保護與生態文明建設提供參考。

1研究范圍與方法

1.1研究區域及數據

大沽河地處山東省膠東半島西部，位于120°07′～120°34′E、36°02′～37°05′N之間。河流發源于招遠市東北部的阜山，自北往南流經招遠市、萊西市、平度市、即墨區、膠州市、城陽區6個市（區），于膠州市營海街道辦事處東營村入膠州灣，干流全長199 km，流域面積6 205 km2。

選取大沽河干流南村水文站1956—2016年逐月降水量、實測徑流量、實測輸沙量資料，采用青島市第三次水資源調查評價成果，對南村水文站控制流域內受水利工程和其他因素的影響而損耗或增加的河川徑流量逐年進行還原計算，求得歷年天然月、年徑流量。該數據系列完整，無缺失。南村水文站位于青島市平度市南村鎮南村，控制流域面積3 724 km2，距入海口62 km，在大沽河流域具有較好的代表性。

1.2研究方法

1.2.1變化趨勢研究

本研究采用線性趨勢法[11]、Mann-Kendall秩相關檢驗法[12]（以下簡稱“MK法”）和改進MK法（包括MMK法[13-14]、PW-MK法[15]、TFPW-MK法[16]）進行趨勢分析。

MK法是世界氣象組織推薦并已被廣泛使用的一種非參數檢驗方法，該方法的優點是不需要待檢序列遵從一定的分布。構建檢驗統計量Z，當Z>0時，存在上升的趨勢；當Z<0時，存在下降的趨勢。當樣本數量n比較大時，Z近似服從標準化正態分布，給定顯著性水平α，則可根據Z與臨界值Z（1-α/2）的比較結果，判定序列趨勢的統計顯著性。

MK法基于序列獨立性假設，而水文序列的自相關性將影響檢驗結果的顯著性水平。改進的MMK法、PW-MK法和TFPW-MK法均對時間序列的相關性進行了處理，減小自相關性對檢驗結果的影響，詳見表1。

1.2.2突變檢驗

本研究采用滑動T檢驗法、累積距平曲線法進行突變檢驗。

滑動T檢驗是通過考察兩個樣本平均值的差異是否顯著來檢驗突變點。對于具有n個樣本的時間序列x，人為設置某一時刻為基準點，基準點前后兩端子序列x1和x2的均值差異如果超過了一定的顯著性水平，可以認為均值發生了質變，有突變發生[17]。該方法的缺點是子序列時段的選擇具有人為性。為避免任意選擇子序列長度造成突變點的漂移，具體使用這一方法時，可以反復變動子序列長度進行比較，提高計算結果的可靠性[18]。

累積距平曲線法是一種直觀判斷趨勢變化的方法，核心是離散數據大于平均值，累積距平值增大，曲線呈現上升趨勢，反之則呈現下降趨勢；根據累積距平值曲線的波動起伏，可以判斷長期的演變趨勢以及變化趨勢發生突變的時間[19]。

1.2.3周期分析

本研究應用Morlet連續復小波分析法[20-21]研究徑流量和輸沙量的周期性變化規律，提取變化的主周期，揭示多時間尺度變化的復雜結構。小波分析是一種窗口大小固定但時間窗和頻率窗都可改變的時頻局部化方法，可以同時揭示時頻細節信息。小波方差隨尺度的變化過程稱小波方差圖，它反映了波動的能量隨尺度的分布情況。通過小波方差圖，可以確定一個水文序列中存在的主要時間尺度，即主周期。

1.2.4影響因素的貢獻率計算

本研究采用雙累積曲線法[22]分析不同時期氣候和下墊面變化對徑流量影響的貢獻率，采用累積量斜率變化率法[23-24]分析氣候變化和人類活動對輸沙量影響的貢獻率。

1）雙累積曲線法。通過突變檢驗識別徑流量時間序列的突變點，以徑流量變異點前作為受人類活動影響較小的基礎期，采用雙累積曲線建立徑流量變異前累積天然徑流量與累積降水量關系式，并用人類活動影響顯著時期（變異期）累積降水量計算得到變異期的模擬累積天然徑流量，認為模擬徑流量與基準期人類活動影響條件近似相同，對比變異期與基礎期的實測徑流量、天然徑流量、模擬徑流量等變化，可分析不同時期降水量和人類活動對徑流量影響的貢獻率。

2）累積量斜率變化率法。設累積輸沙量—年份線性關系式的斜率在變異點前后兩個時期分別為SR1和SR2，則累積輸沙量斜率變化率計算公式[24]為

RS=（SR2-SR1）/SR1×100%（1）

RS為正表示斜率增大，為負表示斜率減少。同理，求得累積降水量斜率變化率PS。

降水量變化對輸沙量變化的貢獻率Cp的計算公式為

CP=PS/RS×100%（2）

人類活動對輸沙量變化的貢獻率CH計算公式為

CH=1-CP（3）

2水沙演化特征

2.1變化趨勢

大沽河南村水文站降水量、徑流量與輸沙量年際變化線性趨勢分析見圖1。由圖1可知，1956—2016年大沽河降水量、實測徑流量、天然徑流量、輸沙量均呈下降趨勢，平均減少速率分別為2.146 3 mm/a、0.118 3億m3/a、0.123 1億m3/a、1.927 3萬t/a；其中徑流量、輸沙量均達到顯著性水平（p＜0.05），而降水量未通過顯著性檢驗。

人類活動對大沽河年輸沙量的影響主要表現在修建水利工程和進行水土流失治理。大沽河流域共12座大中型水庫，除北墅水庫（1974年）、廟埠河水庫（1971年）外，其余10座均建成于1958—1960年。此后，20世紀60年代也是小型水庫、塘壩建設的高峰期。由圖1（b）可知，大沽河輸沙量在1965年以后急劇減少。分布在流域范圍內的水庫、塘壩一方面攔蓄了坡面匯入河道的徑流和泥沙，導致大沽河干流特別是中下游輸沙量大幅減少；另一方面，河道內的徑流被攔蓄后，上游來水失去了沖刷下游河道的動力，從而達到減蝕的效果，這在很大程度上改變了流域的水沙狀況。

采用MK法、MMK法、PW-MK法、TFPW-MK法進行趨勢檢驗（見表2），檢驗統計量Z均小于0，說明4個序列都呈下降趨勢；實測徑流量、輸沙量的MK、MMK統計量與天然徑流量的MK、MMK、PW-MK、TFPW-MK統計量|Z|>1.96，下降趨勢顯著（p＜0.05），其他項目均未通過顯著性檢驗。

2.2突變分析

圖2為大沽河南村水文站徑流量、輸沙量累積距平曲線。由圖2可知，1956—1958年徑流量累積距平值呈下降趨勢，1958—1965年出現明顯的上升趨勢，1965—1976年小幅波動，1976—2016年呈現逐漸減少趨勢。由此可見，實測徑流量、天然徑流量的突變點均為1976年，均發生了由多到少的突變。而輸沙量累積距平曲線大致表現為3個變化階段：1956—1965年呈明顯上升趨勢，1965—1979年緩慢下降，1979—2016年呈線性減少趨勢。綜合分析，輸沙量的突變點為1965年、1979年，均發生了由多到少的突變。

采用滑動T檢驗法對徑流量、輸沙量進行突變分析，步長k分別取2，3，…，20，所得結果見圖3。由圖3可知，k≥7時，實測徑流量通過α=0.01顯著性檢驗，

有一個突變點，突變年份為1976年（圖3取k=7為代表曲線）；k≥16時，天然徑流量通過α=0.01顯著性檢驗，有一個突變點，突變年份為1976年（取k=18作為代表曲線）；k≥8時，輸沙量通過α=0.01顯著性檢驗，有兩個突變點，突變年份分別為1965年、1979年（取k=9作為代表曲線）。降水量未檢出突變點。

綜上可知，累積距平法與滑動T檢驗法的結論基本一致，實測徑流量、天然徑流量突變年份為1976年左右，輸沙量突變年份為1965年左右、1979年左右，均發生了由多到少的突變。徑流量、輸沙量突變年份前、后10 a均值及其變幅見表3。

2.3變化周期

2.3.1天然徑流量周期及震蕩情況分析

圖4（a）為天然徑流量小波系數實部圖，正值表示豐水階段，負值表示枯水階段。1956—2016年間，天然徑流量在年際和年代際上都存在明顯周期變化，包括5～9 a和10～15 a的小尺度信號以及16～30 a的大尺度信號。其中，在5～9 a時間尺度有周期表現，但不明顯，其尺度中心在8 a左右；10～15 a時間尺度上正負相位明顯且較穩定，其尺度中心在12 a左右；16～30 a時間尺度內小波曲線閉合完整，正負相位交替出現，說明其周期表現明顯，其尺度中心在21 a左右。從圖4（a）還可以看出，以上3個尺度的周期變化在整個分析時段前期（1976年以前）較為穩定，1976年以后天然徑流量年際間變幅較小，震蕩周期表現不夠穩定。1976年之前小波系數實部等值線密集，1976年之后變稀，說明進入20世紀80年代以來，大沽河徑流量變幅趨緩，處于低流量期，與前述突變檢驗的結論一致。

圖4（b）為天然徑流量小波方差圖，圖中共有3個峰值，第一峰值為21 a，該時間尺度波動能量最強，是豐枯變化的主周期；其次還存在12 a和8 a兩個峰值，分別對應著第二、第三主周期。上述3個周期共同對天然徑流量變化起作用，3個周期的波動控制著天然徑流量在整個時間域內的變化特征。

2.3.2輸沙量周期及震蕩情況分析

圖5為大沽河南村水文站輸沙量小波系數實部、小波方差圖。由圖5可知，輸沙量呈明顯的周期性變化，包括5～8 a和9～17 a的小尺度信號以及18～30 a的大尺度信號。以上3個尺度的周期變化在整個分析時段前期（1979年以前）表現比較穩定。1979年之前小波系數實部等值線密集，1979年之后變稀，說明進入20世紀80年代以來，大沽河輸沙量變幅趨緩，處于低輸沙量期，與前述突變檢驗的結論一致。由圖5（b）可知，輸沙量變化的第一主周期為12 a，第二、第三主周期分別為21 a和7 a。

3氣候變化與人類活動對大沽河水沙演化的貢獻率

3.1對徑流量變化影響的貢獻率

由上述分析可知，大沽河徑流量在1976年有1個突變點。因此，近似假定1956—1976年間是受人類活動影響相對較小、下墊面變化較小的基準期，建立這一時期的累積降水量—天然徑流量雙累積曲線關系，將1977—2016年累積降水量數據代入雙累積曲線關系式，計算該時段在假設下墊面條件與1956—1976年相近條件下的累積徑流量，從而可得到各年的擬合徑流量R擬合。將1980—1989年、1990—1999年、2000—2009年、2010—2016年各時段的擬合徑流量、天然徑流量、實測徑流量分別取均值，記為R擬合、R天然、R實測。R擬合與R天然、R實測對比關系見表4。

人類活動對地表徑流的影響，可分為水資源開發利用活動的直接影響（主要是從河道內取水、修建水庫閘壩等攔蓄工程）和下墊面漸變累積的間接影響。人類活動對流域徑流量的直接影響可表示為實測徑流量R實測與天然徑流量R天然的差值，即對應表4中的“耗水量（直接影響）”一列；人類活動對徑流量的間接影響可表示為1976年后R天然與R擬合差值，即下墊面的影響，對應表4中的“下墊面影響（間接影響）”一列；人類活動對徑流量的綜合影響為1976年后實測徑流量R實測與R擬合的差值，即表4中“綜合影響”一列。

由表4可見，與氣候變化相比，人類活動是造成大沽河徑流量變化的主要因素，影響貢獻率在76.8%～97.4%之間。人類活動的直接影響和間接影響都促使徑流量減小，人類活動的直接影響在20世紀90年代之前逐漸增大，2000—2016年趨于減小；人類活動的間接影響在1980—2016年期間呈現“增大—減小—增大”的趨勢；人類活動的綜合影響在1990—1999年較大，2000年以后比較穩定。

發生這種情況的主要原因是，20世紀80年代初期膠東半島遭遇大旱，青島市供水頻頻告急，大沽河水資源被大量開采。此后，大沽河流域與周邊區域經濟社會發展加速，河道內取水量不斷增加，徑流量不斷減少，強烈的人類活動促使大沽河基本演變為季節性河流[7]。

3.2對輸沙量變化影響的貢獻率

由上述分析可知，輸沙量有2個突變點，分別為1965年、1979年。據此將研究階段劃分為3個時期：A（1956—1965年）、B（1966—1979年）、C（1980—2016年）。由于大沽河流域在20世紀60年代以前受人類活動的影響很小，因此可將1956—1965年這一時期作為基準期，1966—1979年和1980—2016年兩個時期為計算期。3個不同時期累積輸沙量與年份、累積降水量與年份關系曲線見圖6、圖7。表5給出了累積輸沙量和累積降水量的斜率及其變化情況。

氣候變化對輸沙量變化的貢獻率為累積降水量-年份關系斜率變化率與累積輸沙量-年份關系斜率變化率的比值[23-24]。因此，與基準期（A）相比，計算期（B）的降水量和人類活動對輸沙量減小的貢獻率分別為23.8%和76.2%，計算期（C）的降水量和人類活動對輸沙量減小的貢獻率分別為21.3%和78.7%；與計算期（B）相比，計算期（C）的降水量和人類活動對輸沙量減小的貢獻率分別為1.6%和98.4%。

除水利工程外，水土保持治理措施也是影響大沽河年輸沙量的重要因素。大沽河上游主要為青島平度市與萊西市。新中國成立后，該地區先后開展了農村“四旁”植樹、荒山荒灘造林、山區小流域治理、農田林網建設等，尤其在1962年、1981年兩次林木確權發證后，農民植樹的積極性明顯高漲。1985年，平度市森林覆蓋率為10.7%，20世紀末達到20.78%，2005年為26.55%。1983年，萊西市森林覆蓋率為7.97%，1985年為21.5%，2012年達到34.4%。增加植被覆蓋率，增強下滲能力和蓄水能力，可有效減少流域侵蝕產沙，保持水土，從而減少河道內輸沙量[25]。

4結論

1）1956—2016年大沽河實測徑流量、天然徑流量、輸沙量均呈顯著減少趨勢，平均減少速率分別為0.118 3億m3/a、0.123 1億m3/a、1.927 3萬t/a；降水量呈不顯著減少趨勢，平均減少速率為2.146 3 mm/a。

2）1956—2016年大沽河實測徑流量、天然徑流量有1個突變點（1976年），發生了由多到少的突變，突變年份前、后10 a均值分別減少了3.22億m3、2.56億m3，變幅分別為73%、40%。輸沙量有兩個突變點（1965年、1979年），均發生了由多到少的突變，突變年份前、后10 a均值分別減少了113.00萬t、22.94萬t，變幅分別為84%、91%。

3）大沽河徑流量、輸沙量均呈現明顯的周期性震蕩，天然徑流量第一主周期為21 a，第二、第三主周期分別為12 a和8 a；輸沙量第一主周期為12 a，第二、第三主周期分別為21 a和7 a。

4）人類活動是造成大沽河徑流量、輸沙量減少的主要原因，人類活動引起的下墊面變化對天然徑流量變化的貢獻率為76.8%～97.4%；人類活動對輸沙量變化影響的貢獻率為76.2%～98.4%。

[參考文獻]

［1］ 胡春宏，王延貴，張燕菁，等.中國江河水沙變化趨勢與主要影響因素[J].水科學進展，2010，21（4）：524-532.

［2］ 王赫，朱秀端，林敬蘭，等.汀江流域長汀段水沙變化特征及其影響因素[J].水土保持學報，2022，36（2）：122-129，137.

［3］ ELHAM S，MAJID M，ARTEMI C.Sediment transport mechanisms and selective removal of soil particles under unsteadystate conditions in a sheet erosion system[J].International Journal of Sediment Research，2022，37（2）：151-161.

［4］ 李昌文，黃艷，嚴凌志.變化環境下長江流域超標準洪水災害特點研究[J].人民長江，2022，53（3）：29-43.

［5］ 王亞迪，權全，薛濤濤，等.氣候變化對黃河源區的水文影響分析[J].水文水資源，2018，7（2）：135-143.

［6］ 韓沙沙，談廣鳴，趙連軍，等.黃河入海水沙情勢變化分析[J].中國農村水利水電，2022（2）：1-5，19.

［7］ 姜德娟，王曉利.膠東半島大沽河流域徑流變化特征[J].干旱區研究，2013，30（6）：965-972.

［8］ 盛茂剛，黃修東，左林遠，等.大沽河流域近60年降水量時空變化特征分析[J].水資源與水工程學報，2016，27（2）：65-68.

［9］ 金秉福，宮立新，宋鍵.大沽河泥沙來源的重礦物分析及其環境意義[J].海洋科學，2010，34（10）：71-76.

［10］ 薛偉，徐映雪，羊艷.青島市大沽河地下水水源地應急能力分析[C].//水利部水文局，國際水文計劃中國國家委員會，國際水文科學協會中國國家委員會，等.中國水文科技新發展：2012中國水文學術討論會論文集.南京：河海大學出版社，2012：833-836.

［11］ 趙嘉陽，郭福濤，梁慧玲，等.福建長汀紅壤區1965—2013年氣溫和降水量的變化趨勢[J].福建農林大學學報（自然科學版），2016，45（1）：77-83.

［12］ 張潤潤.香港地區降水趨勢及其演變過程分析[J].河海大學學報（自然科學版），2010，38（5）：505-510.

［13］ AHAMAD I，TANG D S，WANG T F，et al.Precipitation trends over time using Mann-Kendall and Spearmans rho Tests in Swat River Basin，Pakistan[J].Advances in Meteorology，2015（2015）：1-15.

［14］ 白勇.變化環境下半干旱流域徑流演化特征及驅動因素分析[D].呼和浩特：內蒙古農業大學，2018：8-9.

［15］ 章誕武，叢振濤，倪廣恒.基于中國氣象資料的趨勢檢驗方法對比分析[J].水科學進展，2013，24（4）：490-496.

［16］ 許繼軍，楊大文，雷志棟，等.長江流域降水量和徑流量長期變化趨勢檢驗[J].人民長江，2006，37（9）：63-67.

［17］ 馬耘秀，夏露，范毅，等.土地利用和氣候變化對黃土區典型流域水沙變化的影響[J].水土保持學報，2021，35（2）：38-46.

［18］ 魏鳳英.現代氣候統計診斷與預測技術[Ｍ].2版.北京：氣象出版社，2007：63-65.

［19］ 夏偉，周維博，李文溢，等.氣候變化和人類活動對灃河流域徑流量影響的定量評估[J].水資源與水工程學報，2018，29（6）：47-52.

［20］ 魏曉婷，黃生志，黃強，等.定量分解氣候變化與人類活動對季節徑流變異的貢獻率[J].水土保持學報，2019，33（6）：182-189.

［21］ 祁文燕，錢鞠，葛雷，等.湟水干流近60年徑流變化特征分析[J].水資源與水工程學報，2018，29（3）：45-49.

［22］ 師忱，袁士保，史常青，等.灤河流域氣候變化與人類活動對徑流的影響[J].水土保持學報，2018，32（2）：264-269.

［23］ 郭文獻，趙瑞超，付體江，等.烏江流域水沙變化趨勢及驅動力分析[J].人民長江，2021，52（9）：71-78.

［24］ 王隨繼，李玲，顏明.氣候和人類活動對黃河中游區間產流量變化的貢獻率[J].地理研究，2013，32（3）：395-402.

［25］ 彭濤，田慧，秦振雄，等.氣候變化和人類活動對長江徑流泥沙的影響研究[J].泥沙研究，2018，43（6）：54-60.

收稿日期：2022-06-22

基金項目：山東省地礦局海岸帶地質環境保護重點實驗室2021年度開放基金項目（SYS202110）；西北農林科技大學旱區農業水土工程教育部重點實驗室訪問學者基金項目

第一作者：胡萌（1981—），男，山東日照人，高級工程師，博士，主要從事水資源方面的研究。

通信作者：雷炳霄（1984—），男，陜西渭南人，高級工程師，學士，副院長，主要從事水工環地質方面的研究。

E-mail：bingxiao8426@163.com

（責任編輯楊傲秋）