周天祥 陳 晶

(1.紹興文理學院 機械與電氣工程學院，浙江 紹興 312000;2.紹興文理學院 土木工程學院，浙江 紹興 312000)

自然界存在著許多復(fù)雜的動力學過程，這些復(fù)雜的相互作用通過一系列看似嘈雜且不穩(wěn)定的信號表現(xiàn)出來，通?？梢允褂梅蔷€性微分方程來描述這些過程.一些非線性系統(tǒng)的時間序列信號存在著顯著的分形特征[1]，即局部時間序列與整個時間序列完全相同或近似.經(jīng)過許多不同領(lǐng)域?qū)W者的努力，這種現(xiàn)象目前在經(jīng)濟學、生理學、地質(zhì)學以及氣象學等學科的復(fù)雜系統(tǒng)中均有發(fā)現(xiàn)[2-5].大氣中存在很多不規(guī)則的隨機運動(如大氣湍流)，可將其視為在各種尺度的天氣系統(tǒng)相互作用下的復(fù)雜系統(tǒng)，這些相互作用將導致在許多時空尺度上氣象參數(shù)以近似隨機的方式波動.Pielke等人[6]使用非線性微分方程對氣象過程進行建模，但是數(shù)值方法與現(xiàn)實情況仍存在差距.

降水的時空變化對農(nóng)業(yè)生產(chǎn)、水資源管理以及生態(tài)系統(tǒng)有著深刻的影響，同時也是氣候變化的重要參考指標.降雨量依賴于大氣過程的多樣性變化，并由于自然過程的復(fù)雜性在不同的時空尺度上呈現(xiàn)出隨機波動.一些學者基于統(tǒng)計學如Mann-Kendall秩次相關(guān)檢驗、Pettitt檢驗、Buishand檢驗等方法對降水時間序列的變化趨勢進行評價[7-9]，然而這些方法無法捕捉降水時間序列復(fù)雜非線性行為中的潛在特征.隨著非線性動力學的發(fā)展，多重分形級聯(lián)模型逐漸被應(yīng)用于研究天氣與氣候變化.Lovejoy和Mandel-brot[10]發(fā)現(xiàn)降雨區(qū)域的云層邊界投影到地面的形狀為分形圖形，并且降水序列的時空結(jié)構(gòu)服從雙曲線分布，從而確立了分形理論在氣象學的適用性.Telesca等人[11]應(yīng)用去趨勢波動分析方法(DFA)對阿根廷中部地區(qū)近60年的降水時間序列進行研究，發(fā)現(xiàn)了降水時間序列中的標度不變性以及長程相關(guān)性.然而DFA方法適用于描述單個標度指數(shù)，許多復(fù)雜系統(tǒng)中的信號并非表現(xiàn)出簡單的單分形特征，而是在不同的時間尺度上具有不同的分形結(jié)構(gòu)，即多重分形行為.多重分形去趨勢波動分析方法(MF-DFA)能更加精確地描述復(fù)雜系統(tǒng)可能存在的多個分形特征[12].Kalamaras等人[13]應(yīng)用MF-DFA方法對希臘克里特島沿海氣象站的氣溫觀測數(shù)據(jù)進行了研究，發(fā)現(xiàn)該地區(qū)的氣溫時間序列表現(xiàn)出多重分形特征，并且氣溫時間序列為長期正相關(guān).Olsson[14]運用概率密度函數(shù)以及功率譜分析了降水時間序列所服從的冪律分布，選取的不同統(tǒng)計矩階數(shù)隨尺度而變化證明該地區(qū)的降水時間序列具有多重分形特征.新疆北部地區(qū)水資源稀缺，降雨量的變化對該地區(qū)農(nóng)牧業(yè)的生產(chǎn)至關(guān)重要.本文基于新疆北部阿勒泰地區(qū)1961至2015年的逐月降雨量數(shù)據(jù)，對該降水時間序列的分形特征和變化趨勢運用MF-DFA方法進行了詳細分析.

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

本文研究資料來源于額爾齊斯河流域標準氣象站觀測數(shù)據(jù)[15]，研究的三個地點分別位于新疆北部阿爾泰地區(qū)的阿勒泰氣象站(47.73 °N，88.08 °E)、福海氣象站(47.12 °N，87.47 °E)、布爾津氣象站(47.70 °N，86.87 °E)，所在位置已在圖1中標明.該區(qū)域深居歐亞大陸腹地，屬于溫帶大陸性氣候，冬季在強大的蒙古高壓控制下氣候寒冷干燥，年平均降水量在100 mm左右，降水多集中在夏季，冬季降雨量較春秋季節(jié)更多.三個觀測站的降水資料數(shù)據(jù)的時間范圍為1961年1月至2015年12月，觀測數(shù)據(jù)的時間尺度為月.

2 MF-DFA原理

分形是指局部與整體之間具有某種方式相似的形狀，即局部片段與整體之間具有自相似性，目前分形理論已從狹義的幾何對象擴展至廣義分形.Mandelbrot[16]在研究美國棉花價格的變化時發(fā)現(xiàn)價格時間序列的短期波動與長期的變化趨勢存在標度不變性，因而發(fā)現(xiàn)了時間序列中的自相似現(xiàn)象.若分形維隨時間尺度的不同而發(fā)生變化，則該時間序列具備多重分形特性.多重分形去趨勢波動分析方法(MF-DFA)由Kandelhardt等人[12]在消除趨勢波動分析(DFA)方法的基礎(chǔ)上拓展得到，能更加精確地描述復(fù)雜系統(tǒng)的動力學行為，具體計算方法如下：

第一步，首先對降水時間序列xi(i=1，2，…，n)求取累積離差序列：

(1)

第二步，將Yi使用長度為s的Ns=int(N/s)個互不重疊的區(qū)間，對于每個長度為s的區(qū)間進行最小二乘擬合，得到該區(qū)間的局部趨勢多項式y(tǒng)v(s).對Yi上的每個子區(qū)間v計算局部均方誤差F2(v，s)：

圖1 研究區(qū)域所在位置

(2)

第三步，第qth階波動函數(shù)可表示為：

(3)

第四步，對指定的q階值可使用下式確定波動函數(shù)的標度系數(shù)：

Fq(s)∝sH(q)

(4)

式中H(q)為廣義Hurst指數(shù)，其取值范圍為[0，1]，H(q)越接近1，表明時間序列存在的噪聲越少，趨勢越明顯.當H=0.5時，序列變化為標準布朗運動，序列之間不相關(guān)，表明當前降水時間序列的趨勢不會影響未來的情況；當H(q)>0.5時，表明當前降水時間序列與未來增量之間為正相關(guān)；當H(q)<0.5時，表明當前的降水時間序列與未來增量之間為負相關(guān).若H(q)指為常數(shù)時表明時間序列為單分形的情況，當H(q)指數(shù)隨q值的增大而減小時表明時間序列屬于多重分形情況[17].

第五步，計算多重分形標度系數(shù)τ(q)滿足下列關(guān)系式：

τ(q)=qH(q)-1

(5)

第六步，根據(jù)Legendre變換，多重分形奇異譜函數(shù)f(α)滿足下式：

f(α)=qα-τ(q)

(6)

其中α=?τ(q)/?(q)為奇異性指數(shù)，奇異譜的譜寬(Δα=αmax-αmin)描述了多重分形的程度，譜寬越大表明多重分形程度越高，降水時間序列波動越劇烈.若時間序列具有多重分形屬性，其奇異譜將呈現(xiàn)為倒拋物線形[18].Δf=f(αmin)-f(αmax)表示最大概率子集與最小概率子集的分形差異，Δf可以反映時間序列的變化趨勢，Δf<0表明未來時間序列數(shù)據(jù)的向低數(shù)值方向發(fā)展的概率較大[19].

3 實驗結(jié)果分析

本文選取新疆阿勒泰地區(qū)的阿勒泰、福海、布爾津三個氣象觀測站的逐月觀測數(shù)據(jù)，每個降水觀測數(shù)據(jù)的時間序列長度都是660.由圖2可知，三個地區(qū)的歷年逐月降雨量基本都不足100 mm，降水時間序列全都表現(xiàn)出非線性、非平穩(wěn)性的特點.福海地區(qū)降雨量最少，且在1984年出現(xiàn)一次大降水天氣，降水量超過了100 mm.阿勒泰地區(qū)降雨量趨勢受季節(jié)變化影響明顯，夏季降雨量要明顯高于其他季節(jié)，因此應(yīng)先將原始的降水時間序列轉(zhuǎn)化為隨機游走結(jié)構(gòu)[20]以剔除季節(jié)更替帶來的影響.通過對三個氣象站的降水時間序列xi計算累積離差(式2)可得到隨機游走的降水時間序列Yi，如圖3所示，總體上降水累積離差序列在20世紀60年代至90年代呈下降趨勢，而在60年代末和80年代出現(xiàn)了兩個小幅震蕩的時期，90年代后又呈增加趨勢.

(a)阿勒泰 (b)福海 (c)布爾津

(a)阿勒泰 (b)福海 (c)布爾津

對三個地區(qū)降水隨機游走時間序列Yi按照時間尺度s進行分割，為防止多項式過擬合降水時間序列從而導致剩余波動接近于0，時間尺度s應(yīng)滿足不小于多項式擬合階數(shù)，同時最大時間尺度應(yīng)是時間序列劃分區(qū)間的個數(shù)大于10個[21].統(tǒng)計矩階數(shù)q的取值范圍對多重分形譜的計算十分關(guān)鍵，q值過大會導致計算時間呈指數(shù)增大，q值過小時多重分形譜的波動會很大[22].由此本文選取Yi的子區(qū)間長度s=(6∶66)，按照經(jīng)驗取q=(-10∶10)，利用式(3)對不同取值的統(tǒng)計矩階數(shù)q計算波動函數(shù)Fq(s).波動函數(shù)Fq(s)關(guān)于時間尺度s在特定統(tǒng)計矩階數(shù)q=(-10，0，10)下的雙對數(shù)圖像如圖4所示，所有q值情況下的波動函數(shù)Fq(s)關(guān)于時間尺度s都近似為線性關(guān)系，這表明三個地區(qū)的降水時間序列都符合冪律分布，同時也證明了降水時間序列的內(nèi)部波動演化并非隨機行為，而是由遠距離的某些特征有序引發(fā).另外，從圖4可以明顯看出，回歸線的斜率明顯隨q值的增大而減小，由式(4)可知該斜率即為廣義Hurst指數(shù)，圖4中的擬合方程以及可決系數(shù)見表1.q值描述了局部波動均方根占整體的比重[21]，當q值為正時，廣義Hurst指數(shù)表示大波動的標度行為，q值為負時，廣義Hurst指數(shù)表示小波動的標度行為，對于較小的時間尺度s受降水時間序列波動幅度影響較大，而在較大的時間尺度上降水時間序列的波動幅度將被削弱[22].由圖2(b)可以看出，福海地區(qū)的降水時間序列出現(xiàn)過一次大降水天氣情況，同時其累積離差分布也存在多個波動趨勢，降水時間序列受大的波動控制，因此其廣義Hurst指數(shù)變化范圍最大(ΔH(q)=0.57).

(a)阿勒泰 (b)福海 (c)布爾津

廣義Hurst指數(shù)可用于判斷降水時間序列的分形特征以及長程相關(guān)性.如圖5所示，三個地區(qū)的Hurst指數(shù)關(guān)于q為非線性關(guān)系，說明三個地區(qū)的降水時間序列都存在多重分形特征.當q值為2時，廣義Hurst指數(shù)與經(jīng)典的Hurst指數(shù)相同，此時阿勒泰、福海、布爾津地區(qū)降水時間序列的廣義Hurst指數(shù)分別為0.56、0.50、0.52，因此當q∈[-10，2]時，H(q)>0.5，三個地區(qū)的降水時間序列都具有長期正相關(guān)性.由圖5可知，阿勒泰與布爾津的降水時間序列的長程相關(guān)性明顯強于其他福海地區(qū)，表明未來的該地區(qū)降水趨勢將與當前保持一致.

表1 降水時間序列廣義Hurst指數(shù)的計算

圖5 降水時間序列的廣義Hurst指數(shù)

根據(jù)(式11)計算三個地區(qū)降水時間序列的多重分形譜如圖6所示，其多重分形譜均呈倒拋物線形，分形譜的最大值對應(yīng)為q=0時.由表2可知，福海地區(qū)的譜寬最大，Δα=0.96，因此該地區(qū)降水時間序列的多重分形程度要強于阿勒泰以及布爾津地區(qū)，并且其降水時間序列波動幅度最劇烈.阿勒泰氣象站和布爾津氣象站的譜寬比較接近，因此其多重分形程度也相似.阿勒泰和福海的降水時間序列的譜高度Δf<0，表明這兩個地區(qū)未來降雨量可能有降低的趨勢.布爾津地區(qū)降水時間序列的譜高度Δf>0，因此未來該地區(qū)降雨量增加的概率較大.

表2 降水時間序列的多重分形參數(shù)

(a)阿勒泰 (b)福海 (c)布爾津

4 結(jié)論

本文通過對新疆北部阿勒泰地區(qū)的阿勒泰、福海、布爾津三個氣象站的1961年至2015年的逐月降雨量時間序列的變化進行多重分形檢驗，結(jié)論如下：

(1)根據(jù)MF-DFA方法的檢驗結(jié)果，三個地區(qū)的降水資料時間序列都表現(xiàn)出顯著的多重分形特征，符合冪律相關(guān)特性；

(2)計算廣義Hurst指數(shù)發(fā)現(xiàn)降水時間序列在q<2時存在明顯的長期正相關(guān)，表明未來的降水趨勢具備可預(yù)測性；

(3)福海地區(qū)的降水時間序列相比另外兩個地區(qū)的多重分形譜寬更大、波動更強烈，表現(xiàn)出更強的多重分形程度；

(4)對于長期的降水時間序列使用MF-DFA方法，可以獲取其統(tǒng)計學的表征以便更好地理解其動力學機制，為未來當?shù)亟邓淖兓厔萏峁﹨⒖?