王俊杰, 拾 兵, 巴彥斌

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

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

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

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

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

2 研究方法

趨勢(shì)性、周期性和突變性是水文時(shí)間序列的重要特征[13]，本文采用5 a滑動(dòng)平均法、Mann-Kendall突變檢驗(yàn)法、滑動(dòng)T檢驗(yàn)等方法對(duì)水沙序列的趨勢(shì)性和變異點(diǎn)進(jìn)行分析，并通過小波變換、經(jīng)驗(yàn)?zāi)B(tài)分解等方法對(duì)水沙通量的周期性進(jìn)行研究，利用交叉小波分析和小波相干譜對(duì)徑流量和輸沙量的共振周期與相干性進(jìn)行研究，其中Mann-Kendall突變檢驗(yàn)法、滑動(dòng)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)信號(hào)的分解方法，它利用信號(hào)內(nèi)部的特征尺度變化作頻率與能量的解析，將非平穩(wěn)、非線性信號(hào)分解為若干本征模態(tài)函數(shù)(intrinsic mode function,IMF)和一個(gè)殘余分量。對(duì)于給定的原始信號(hào)X(t)，其分解過程為：

(1)

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

(2)

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

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

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

(3)

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

(4)

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

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

(5)

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

3 結(jié)果與分析

3.1 水沙通量趨勢(shì)及變異分析

為分析黃河入海水沙通量的趨勢(shì)，分別對(duì)1950—2018年利津站水沙時(shí)間序列做趨勢(shì)分析和五年滑動(dòng)平均處理。近70 a黃河入海年徑流量呈現(xiàn)波動(dòng)下降的趨勢(shì)見圖2A。在1950—1970年，徑流量平均值較高，處在波動(dòng)的偏高期，1970年之后整體下降較快；近70 a下降的趨向率為64.3億m3/10 a，年均徑流量為293.7億m3。輸沙量在近70 a也呈現(xiàn)出與徑流量相似的趨勢(shì)特征見圖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，在人類活動(dòng)(水庫(kù)修建、水土保持措施、干流引水等)和氣候變化雙重因素的作用下[16]，黃河入海水沙通量發(fā)生了深刻變化。1970年之前，黃河流域大型水利樞紐工程較少，河流連通性順暢，且流域降雨充沛，入海水沙量主要受氣候變化影響；70年代以后，黃河干支流陸續(xù)修建了劉家峽、龍羊峽、青銅峽、小浪底等水利樞紐，黃河流域的徑流逐漸被干支流的水庫(kù)控制，其中，黃河干流中、上游水庫(kù)除三門峽外共有7座，總庫(kù)容高達(dá)314.12億m3[11]，改變了黃河天然的徑流狀態(tài)，河流挾沙能力下降，1986年之后的年平均入海泥沙僅為1986年之前的23.2%。

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

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

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

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

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

3.2 水沙通量多時(shí)間尺度特征

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

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

圖5 徑流量EMD分解

圖6 輸沙量EMD分解

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

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

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

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

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

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

3.3 水沙通量的共振周期

對(duì)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)系；由徑流和輸沙的小波相干譜可知，黃河入海徑流和輸沙在整個(gè)時(shí)頻空間低能量區(qū)內(nèi)存在顯著的共振周期，且低能量區(qū)的顯著相關(guān)性遠(yuǎn)大于高能量區(qū)，其中，8～16 a的時(shí)間尺度在整個(gè)時(shí)間段內(nèi)均通過95%的顯著性檢驗(yàn)(兩端時(shí)間受邊界影響)，且二者以正相位為主，說(shuō)明黃河入海徑流和輸沙演變特征具有極強(qiáng)的一致性。

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

4 結(jié) 論

(1) 近70 a黃河入海水沙呈顯著下降趨勢(shì)，黃河入海徑流量在1985年發(fā)生變異，輸沙量在1996年發(fā)生變異，且在1970年之后，黃河入海水沙通量的各種尺度信號(hào)出現(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的時(shí)間尺度上具有共同的振蕩周期，顯著性區(qū)域集中在1960—1972年，表明黃河入海水沙之間的相關(guān)性在該時(shí)頻空間下最為顯著。小波相干譜分析表明，輸沙和徑流以正相位關(guān)系為主，且黃河入海水沙之間的低能量區(qū)的顯著性遠(yuǎn)大于高能量區(qū)，反映出黃河入海徑流和輸沙的演變特征具有一致性。