王俊杰, 拾 兵, 巴彥斌

(中國海洋大學(xué) 工程學(xué)院, 山東 青島 266100)

河流入海水沙通量對河口地貌塑造、岸灘侵蝕以及近海岸泥沙輸移有著重要作用[1]；同時，在河流水動力的作用下會使得泥沙顆粒攜帶大量污染物遷移到河口及近海岸，嚴(yán)重影響河口生態(tài)環(huán)境[2]；再者，河口及近海岸區(qū)域存在的懸浮泥沙顆粒會與水體微生物以及其他有機(jī)物相互作用形成懸浮物，懸浮物的大量存在會影響水體透明度，進(jìn)一步影響水體生物的光合作用和初級生產(chǎn)力[3]。因而關(guān)于入海水沙通量的研究是陸海相互作用研究的重要科學(xué)問題之一。

目前，關(guān)于黃河水沙過程的研究已引起了廣泛關(guān)注，但研究區(qū)域多集中在中上游段，并且對水沙過程的特征分析不夠全面。如張金萍等[4]僅研究了黃河潼關(guān)站水沙過程的趨勢性、周期性，并未對其變異性進(jìn)行討論，且方法較為單一。歐陽潮波等[5]對黃河河龍區(qū)間的水沙長期趨勢、年際變異、以及水沙過程對人類活動的響應(yīng)進(jìn)行了討論，對水沙的周期性并未涉及。姚文藝等[6]利用“水文法”和“水保法”對黃河上中游水沙演變的驅(qū)動成因和未來趨勢進(jìn)行了研究，認(rèn)為未來黃河水沙量仍處于下降趨勢。另外，其他眾多研究成果[7-10]表明黃河入海水沙通量呈顯著下降的趨勢，且人類活動(水庫建設(shè)、引水引沙和水土保持措施等)的加劇以及氣候變化的影響是水沙通量減少的主要因素。那么，在多種耦合因素的影響下，最近70 a黃河入海水沙呈現(xiàn)怎樣的多尺度特征？徑流和輸沙是否存在相似的演變規(guī)律？以及水沙之間是否具有協(xié)同影響和共振周期？為了明晰以上問題，本文以黃河利津站1950—2018年水沙資料為基礎(chǔ)，應(yīng)用滑動平均法對水沙趨勢性進(jìn)行討論，以Mann-Kendall突變檢驗(yàn)法、滑動T檢驗(yàn)等對水沙序列進(jìn)行年際變異性分析，同時利用經(jīng)驗(yàn)?zāi)B(tài)分解和Morlet小波分析對水沙過程的多尺度周期進(jìn)行探討，另外利用交叉小波分析和小波相干性對徑流和輸沙的共振周期進(jìn)行研究。旨在為進(jìn)一步加強(qiáng)對黃河口及其鄰近海域的生態(tài)治理、三角洲的侵蝕保護(hù)等提供理論依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)來源

黃河流域處于東亞海陸季風(fēng)區(qū)的北部，是中國的第二大河，以水少沙多而聞名于世，具有“水沙異源”的顯著特點(diǎn)[11]。利津站是黃河流域的入?？刂普?在利津水文站至入?？谥g的泥沙淤積量無法確定時,可將利津站實(shí)測懸移質(zhì)輸沙量近似視為黃河入海泥沙通量[12]。本文所探討的黃河入海水沙通量數(shù)據(jù)為1950—2018年利津水文站實(shí)測資料，其主要來源分為三部分：一是《全國主要河流水文特征統(tǒng)計(jì)》，時間跨度為1950—1979年；二是黃河利津水文站，時間跨度為1980—1999年；三是《中國河流泥沙公報(bào)》，時間跨度為2000—2018年。資料均通過水文的“三性”審查，研究區(qū)利津水文站地理位置如圖1所示。

2 研究方法

趨勢性、周期性和突變性是水文時間序列的重要特征[13]，本文采用5 a滑動平均法、Mann-Kendall突變檢驗(yàn)法、滑動T檢驗(yàn)等方法對水沙序列的趨勢性和變異點(diǎn)進(jìn)行分析，并通過小波變換、經(jīng)驗(yàn)?zāi)B(tài)分解等方法對水沙通量的周期性進(jìn)行研究，利用交叉小波分析和小波相干譜對徑流量和輸沙量的共振周期與相干性進(jìn)行研究，其中Mann-Kendall突變檢驗(yàn)法、滑動T檢驗(yàn)、小波分析等方法為水文分析常用方法，因此不再贅述，經(jīng)驗(yàn)?zāi)B(tài)分解方法和交叉小波分析方法如下：

圖1 研究區(qū)及水文站位置

(1) EMD方法。經(jīng)驗(yàn)?zāi)B(tài)分解是Huang等人[14]在1998年提出的一種自適應(yīng)信號的分解方法，它利用信號內(nèi)部的特征尺度變化作頻率與能量的解析，將非平穩(wěn)、非線性信號分解為若干本征模態(tài)函數(shù)(intrinsic mode function,IMF)和一個殘余分量。對于給定的原始信號X(t)，其分解過程為：

(1)

② 計(jì)算h1(t)=X(t)-m1(t)

(2)

③ 判斷h1是否是IMF，若不是，則繼續(xù)重復(fù)① 、②過程，直到均值包絡(luò)線趨于零，此時，記h1為c1，c1即為篩選出的第一個IMF，它表示原始信號的最高頻分量。

④ 將c1從原始信號中分離，得到去除一個高頻分量的新信號r1(t)，即：

r1(t)=X(t)-c1(t)

(3)

r1(t)也稱為殘余信號，再對新信號重復(fù)以上步驟，直到第n個殘余信號變?yōu)閱握{(diào)函數(shù)，即可終止分解，最終得到一組IMF和一個殘差信號rn(t)。原始信號X(t)可以重構(gòu)為：

(4)

式中：rn(t)為殘差項(xiàng)，代表信號的平均走勢。

小波相干性(wavelet coherence，WTC)用來反映兩個時間序列在時頻空間的相干程度，其定義為：

(5)

小波相干性顯著性檢驗(yàn)采用以紅噪聲為標(biāo)準(zhǔn)譜的Monte Carlo方法。

3 結(jié)果與分析

3.1 水沙通量趨勢及變異分析

為分析黃河入海水沙通量的趨勢，分別對1950—2018年利津站水沙時間序列做趨勢分析和五年滑動平均處理。近70 a黃河入海年徑流量呈現(xiàn)波動下降的趨勢見圖2A。在1950—1970年，徑流量平均值較高，處在波動的偏高期，1970年之后整體下降較快；近70 a下降的趨向率為64.3億m3/10 a，年均徑流量為293.7億m3。輸沙量在近70 a也呈現(xiàn)出與徑流量相似的趨勢特征見圖2B。近70 a泥沙下降的趨向率為2.14億t/10 a，年輸沙量平均值為6.62億t，兩者的相關(guān)性表明，在α=0.01的置信度水平上呈顯著正相關(guān)，Pearson相關(guān)系數(shù)為0.853。

近70 a，在人類活動(水庫修建、水土保持措施、干流引水等)和氣候變化雙重因素的作用下[16]，黃河入海水沙通量發(fā)生了深刻變化。1970年之前，黃河流域大型水利樞紐工程較少，河流連通性順暢，且流域降雨充沛，入海水沙量主要受氣候變化影響；70年代以后，黃河干支流陸續(xù)修建了劉家峽、龍羊峽、青銅峽、小浪底等水利樞紐，黃河流域的徑流逐漸被干支流的水庫控制，其中，黃河干流中、上游水庫除三門峽外共有7座，總庫容高達(dá)314.12億m3[11]，改變了黃河天然的徑流狀態(tài)，河流挾沙能力下降，1986年之后的年平均入海泥沙僅為1986年之前的23.2%。

在分析1950—2018年黃河入海水沙趨勢的基礎(chǔ)上，通過M-K非參數(shù)檢驗(yàn)法、滑動T檢驗(yàn)、累積距平法對徑流和輸沙的突變年份進(jìn)行確定。在入海徑流量M-K統(tǒng)計(jì)量曲線中，在0.05置信度水平內(nèi)，UF和UB曲線在1985年出現(xiàn)交點(diǎn)見圖3A；入海泥沙量則在1995年，1996年兩個年份出現(xiàn)交點(diǎn)見圖3B。因此，對應(yīng)的時間節(jié)點(diǎn)可能為水沙時間序列的突變年份。

圖2 黃河入海徑流量、輸沙量變化趨勢

圖3 黃河入海徑流量、輸沙量M-K檢驗(yàn)

為了進(jìn)一步闡明近70 a入海輸沙量和徑流量的突變情況，結(jié)合水沙序列的累積距平曲線(圖4)，發(fā)現(xiàn)1950—2018年水沙時間序列的累積距平曲線大致呈倒“V”型，同時其拐點(diǎn)出現(xiàn)在1985年，說明在1950—1985年，入海水沙量呈現(xiàn)增長的趨勢，在1985年之后開始出現(xiàn)下降，整體上經(jīng)歷了由“豐”到“枯”的變化；再進(jìn)一步結(jié)合滑動T檢驗(yàn)發(fā)現(xiàn)，在步長為3，0.01顯著水平下，徑流量UF和UB曲線相交的1985年通過滑動T檢驗(yàn)；對于輸沙量，只有1996年通過了步長為4，顯著水平為0.01的滑動T檢驗(yàn)，因此，在結(jié)合多種方法的基礎(chǔ)上，認(rèn)為入海徑流量在1985年發(fā)生變異，入海輸沙量在1996年發(fā)生變異，盡管徑流量和輸沙量的時間序列具有較強(qiáng)的相關(guān)性(相關(guān)系數(shù)為0.853)，但兩者的變異年份仍有一定時差，這與黃河流域“水沙異源”特性以及人類活動干擾等因素有關(guān)[12]。

圖4 徑流和輸沙的累積距平曲線

3.2 水沙通量多時間尺度特征

本文采用經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)與Morlet小波分析的方法對1950—2018年黃河入海水沙時間序列進(jìn)行多時間尺度特征分析。經(jīng)驗(yàn)?zāi)B(tài)分解是一種數(shù)據(jù)驅(qū)動算法，具有較好的靈活性和自適應(yīng)性，能夠根據(jù)數(shù)據(jù)固有特征進(jìn)行分解[17]，圖5—6分別為徑流量和輸沙量的經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果，可以看出1950—2018年黃河入海徑流量分成了3個本征模態(tài)函數(shù)(IMF)和1個趨勢項(xiàng)(Res)，輸沙量被分解成了4個IMF和1個趨勢項(xiàng)，分別對徑流量和輸沙量的IMF和趨勢項(xiàng)進(jìn)行重構(gòu)，發(fā)現(xiàn)與原始序列完全重合?？梢哉J(rèn)定EMD方法分解結(jié)果正確可靠，同時通過對徑流量和輸沙量的趨勢項(xiàng)(Res)的大致走向分析，發(fā)現(xiàn)徑流量在1950—2000年處于減少的趨勢，在2000年左右以后，出現(xiàn)了小幅度上升態(tài)勢，而輸沙量一直處于減少的趨勢，這與上文的趨勢分析中的結(jié)果基本一致。

為了分析徑流和輸沙的不同時間尺度特征，對徑流和輸沙經(jīng)過經(jīng)驗(yàn)?zāi)B(tài)分解的各個分量利用FFT求周期的方法進(jìn)行分析[18]，結(jié)果見表1。

圖5 徑流量EMD分解

圖6 輸沙量EMD分解

表1 周期識別結(jié)果

項(xiàng)目IMF1IMF2IMF3IMF4徑流3～59～1119～21—輸沙3～45～69～1119～21

由表1可知，徑流和輸沙在不同時間尺度上的周期基本一致，其中輸沙量的IMF1和IMF2的周期尺度相近，這可能由經(jīng)驗(yàn)?zāi)B(tài)分解的局部篩選性質(zhì)制約而造成“模態(tài)混疊”現(xiàn)象有關(guān)[19]，但不影響對水沙時間序列周期的識別。黃河入海徑流量具有3個模態(tài)，第一階模態(tài)IMF1頻率最高，周期最小，IMF1和IMF2的分別具有3～5 a，9～11 a的尺度周期，反映了徑流量的年際尺度振蕩，而IMF3的周期為19～21 a，是入海徑流量年代際尺度周期；黃河入海泥沙具有4個模態(tài)，其中，IMF1是頻率最高，周期最小的模態(tài)，與徑流的IMF1對應(yīng)，具有3～4 a的尺度周期，IMF2，IMF3周期分別為5～6 a，9～11 a，以上3個模態(tài)均為輸沙量的年際尺度周期，IMF4周期為19～21 a，是入海輸沙量的年代際周期，黃河入海水沙呈現(xiàn)出的多尺度特征，表明黃河流域受到人類活動、氣候、下墊面變化等因素的綜合作用，具有復(fù)雜的演變過程。

在采用經(jīng)驗(yàn)?zāi)B(tài)分解對黃河入海水沙進(jìn)行多尺度特征分析的基礎(chǔ)上，亦利用Morlet小波分析進(jìn)行對比，圖7分別為1950—2018年黃河徑流量和輸沙量的小波系數(shù)等值線圖，可以反映不同時間段下的周期振蕩強(qiáng)弱程度，以及水沙時間序列的豐、枯交替的變化現(xiàn)象，其中，實(shí)線為豐水豐沙階段，虛線為枯水枯沙階段。黃河入海徑流量分別存在3～5 a，9～11 a以及20～22 a的尺度周期，與經(jīng)驗(yàn)?zāi)B(tài)分解的結(jié)果較為吻合，并且入海徑流量在年代際尺度周期上大致經(jīng)歷了“豐—枯—豐—枯—豐”5個階段，未來一段時間內(nèi)仍為豐水階段，且20～22 a的尺度周期具有全局性，存在徑流的整個階段，3～5 a，9～11 a的年際尺度周期在1985年之前振蕩較強(qiáng)，在1985年之后逐漸減弱，這與上文M-K檢驗(yàn)、滑動T檢驗(yàn)認(rèn)為徑流在1985年發(fā)生變異的結(jié)論具有一致性；同時，通過分析黃河入海泥沙量的小波系數(shù)等值線圖發(fā)現(xiàn)，1950—2018年入海泥沙量大致存在3個不同的時間尺度周期，分別為2～5 a，9～11 a，22～24 a，輸沙量小波分析的結(jié)論與經(jīng)驗(yàn)?zāi)B(tài)分解的結(jié)果具有一致性，輸沙量在22～24 a的年代際尺度周期上大致經(jīng)歷了“豐—枯—豐—枯”4個階段，并且未來幾年仍為枯沙階段，輸沙量的年代際周期在1970年之后振蕩逐漸減弱，9～11 a尺度周期在整個時間段內(nèi)較為穩(wěn)定，具有全局性。

圖7 徑流量、輸沙量小波系數(shù)等值線

綜合經(jīng)驗(yàn)?zāi)B(tài)分解與小波分析，發(fā)現(xiàn)黃河入海水沙的長期演變過程具有顯著的年際以及年代際尺度周期變化特征，同時，入海徑流量和輸沙量的趨勢性以及周期尺度特征基本吻合。并且，在1985年以后，黃河入海水沙的各種尺度周期振蕩強(qiáng)度逐漸減弱，這與黃河流域大規(guī)模人類活動有關(guān)，黃河入海水沙進(jìn)入了新的階段。

3.3 水沙通量的共振周期

對1950—2018年黃河入海水沙序列進(jìn)行交叉小波分析，圖8徑流量和輸沙量的交叉小波譜和小波相干譜；其中，箭頭方向表示兩者之間的相位關(guān)系，由左向右表示為同相位關(guān)系，由右向左表示反相位，黑色粗實(shí)線表示兩者之間達(dá)到了95%的紅噪聲檢驗(yàn)，黑色的細(xì)實(shí)線為小波影響錐線(COI)，該曲線以外區(qū)域由于受到邊緣效應(yīng)而不予考慮[20]。由徑流和輸沙的交叉小波譜可以發(fā)現(xiàn)，徑流和輸沙共同的高能量周期為1～6 a，6～15 a,并且在1960—1972年1～6 a的周期尺度通過了95%顯著性檢驗(yàn)，且二者呈正相位關(guān)系；由徑流和輸沙的小波相干譜可知，黃河入海徑流和輸沙在整個時頻空間低能量區(qū)內(nèi)存在顯著的共振周期，且低能量區(qū)的顯著相關(guān)性遠(yuǎn)大于高能量區(qū)，其中，8～16 a的時間尺度在整個時間段內(nèi)均通過95%的顯著性檢驗(yàn)(兩端時間受邊界影響)，且二者以正相位為主，說明黃河入海徑流和輸沙演變特征具有極強(qiáng)的一致性。

圖8 徑流和輸沙的交叉小波普與小波相干譜

4 結(jié) 論

(1) 近70 a黃河入海水沙呈顯著下降趨勢，黃河入海徑流量在1985年發(fā)生變異，輸沙量在1996年發(fā)生變異，且在1970年之后，黃河入海水沙通量的各種尺度信號出現(xiàn)逐漸減弱的特征。

(2) 1950—2018年黃河入海水沙具有顯著的周期變化規(guī)律，其中，徑流量具有3～5 a，9～11 a，19～21 a年的顯著周期，輸沙量具有3～4 a，5～6 a，9～11 a，19～21 a的顯著周期變化規(guī)律。

(3) 交叉小波分析表明1950—2018年黃河入海水沙在1～6 a的時間尺度上具有共同的振蕩周期，顯著性區(qū)域集中在1960—1972年，表明黃河入海水沙之間的相關(guān)性在該時頻空間下最為顯著。小波相干譜分析表明，輸沙和徑流以正相位關(guān)系為主，且黃河入海水沙之間的低能量區(qū)的顯著性遠(yuǎn)大于高能量區(qū)，反映出黃河入海徑流和輸沙的演變特征具有一致性。