趙曉東, 張宏建, 周洪亮

1 工業(yè)控制技術(shù)國家重點(diǎn)實(shí)驗(yàn)室, 浙江大學(xué), 杭州 310007 2 中國計(jì)量學(xué)院計(jì)量測(cè)試工程學(xué)院, 杭州 310018

基于混沌理論的河流藻類生長(zhǎng)特性分析
——以德國易北河為例

趙曉東1,2,*, 張宏建1, 周洪亮1

1 工業(yè)控制技術(shù)國家重點(diǎn)實(shí)驗(yàn)室, 浙江大學(xué), 杭州 310007 2 中國計(jì)量學(xué)院計(jì)量測(cè)試工程學(xué)院, 杭州 310018

針對(duì)藻類生長(zhǎng)具有高度非線性特征和實(shí)際采樣樣本間隔稀疏的問題，采用了混沌理論對(duì)采樣序列的混沌特征量進(jìn)行估計(jì)。采用C-C方法估計(jì)時(shí)間序列的嵌入維和延遲時(shí)間，采用G-P算法對(duì)關(guān)聯(lián)維進(jìn)行估計(jì)，并采用小數(shù)據(jù)量法估計(jì)最大Lyapunov指數(shù)，最終可實(shí)現(xiàn)對(duì)最長(zhǎng)可預(yù)測(cè)時(shí)長(zhǎng)的估計(jì)。以易北河為例，對(duì)易北河水體葉綠素a 1997年至2001年間各年的觀測(cè)序列進(jìn)行了混沌分析，分析結(jié)果表明，各年的葉綠素a觀測(cè)序列均具有低維混沌特性，關(guān)聯(lián)維D=2.75—4.02，各年的葉綠素a序列的最長(zhǎng)預(yù)測(cè)時(shí)間變化范圍為8.01—18.94 d，平均為13.98d(約2周)。采用同樣方法對(duì)5a易北河連續(xù)日徑流量時(shí)間序列分析表明，該徑流量也具有低維混沌特性(最大Lyapunov指數(shù)λ1=0.0125)，徑流量的最長(zhǎng)預(yù)測(cè)時(shí)間估計(jì)約為80 d。氣候因素的混沌特性對(duì)藻類生長(zhǎng)表現(xiàn)出的混沌特征的影響可能要大于徑流量等水文因素的影響。

藻類生長(zhǎng); 混沌; 關(guān)聯(lián)維數(shù); Lyapunov指數(shù)

水體藻種濃度指標(biāo)實(shí)際上是多藻種濃度指標(biāo)的綜合體現(xiàn)，受限于藻種種群間生長(zhǎng)與競(jìng)爭(zhēng)特性。Huisman等[1]給出了多藻種在不同營養(yǎng)條件下的數(shù)值模擬，藻種生長(zhǎng)顯示出具有混沌特性，藻種顯示出平衡條件下的多樣性。文中只針對(duì)恒定營養(yǎng)鹽輸入進(jìn)行了數(shù)值模擬，并未給出試驗(yàn)驗(yàn)證結(jié)果。自1978年，Connell等[2]提出非平衡營養(yǎng)條件下的 “中等干擾假說”后，對(duì)水生環(huán)境非平衡營養(yǎng)輸入條件下引發(fā)的藻種多樣性進(jìn)行了大量的數(shù)值模擬和實(shí)驗(yàn)研究[2- 8]。國內(nèi)研究者在對(duì)兩種常見水華藻類進(jìn)行 “脈沖”營養(yǎng)鹽輸入條件下的實(shí)驗(yàn)室培養(yǎng)的研究過程中，發(fā)現(xiàn)具有相同初始生長(zhǎng)條件的藻種在非穩(wěn)態(tài)營養(yǎng)輸入條件下顯示出截然不同的競(jìng)爭(zhēng)特性，在生物量上表現(xiàn)出交替占優(yōu)的特征[9]。對(duì)于實(shí)際河流藻類生長(zhǎng)觀測(cè)序列而言，河流下游某觀測(cè)點(diǎn)獲得的藻類生長(zhǎng)數(shù)據(jù)與上游藻類生長(zhǎng)具有相關(guān)性。影響其相關(guān)性的因素除沿途期間的氣象因素之外，河流徑流量也是主要影響因素之一[10]。可見，河流中藻類這種生長(zhǎng)不確定特征除營養(yǎng)鹽輸入因素影響之外，還包括水文和氣象條件等因素的影響。在影響藻類生長(zhǎng)的各種因子中，早在1987年，Tsonis和Elsner就證明了氣候系統(tǒng)是混沌系統(tǒng)[11]，對(duì)于河流徑流量以及降雨量等水文數(shù)據(jù)的研究也表明具有混沌特征[12- 15]。

在對(duì)地表水體中藻類生長(zhǎng)機(jī)理進(jìn)行描述時(shí)，水體網(wǎng)格的劃分細(xì)密程度，各“藕聯(lián)”數(shù)據(jù)的測(cè)量精度都會(huì)對(duì)最終藻類預(yù)測(cè)產(chǎn)生影響，而當(dāng)該系統(tǒng)具有混沌特性時(shí)，上述計(jì)算或測(cè)量誤差都會(huì)對(duì)藻類長(zhǎng)期預(yù)測(cè)結(jié)果產(chǎn)生極大影響，最終導(dǎo)致不可預(yù)測(cè)。另外，由于藻類指標(biāo)實(shí)際采樣樣本間隔稀疏(理論上應(yīng)無窮小)以及藻類生長(zhǎng)本身具有高度非線性特征，有必要對(duì)藻類實(shí)測(cè)時(shí)間序列是否具有混沌特性進(jìn)行分析。目前，盡管各種預(yù)測(cè)方法都針對(duì)特定湖庫和河流能給出了較為成功的應(yīng)用，但對(duì)于究竟能最大預(yù)測(cè)多長(zhǎng)時(shí)間國內(nèi)外并未有文獻(xiàn)給出這方面的報(bào)道。因此，對(duì)于河流藻類預(yù)測(cè)而言，有必要對(duì)其觀測(cè)序列是否混沌進(jìn)行研究，實(shí)現(xiàn)對(duì)其可預(yù)測(cè)時(shí)長(zhǎng)的估計(jì)。

本文對(duì)混沌理論和所采用的經(jīng)典的混沌分析算法進(jìn)行了說明，并以易北河為例，對(duì)葉綠素a觀測(cè)序列和日徑流量的時(shí)間序列進(jìn)行了混沌特征量估計(jì)，以此判斷是否該時(shí)間序列具有混沌特征。并在此基礎(chǔ)上，對(duì)該時(shí)間序列的最大可預(yù)測(cè)時(shí)間進(jìn)行了估計(jì)，并指出了可能成為影響藻類生長(zhǎng)混沌特性的主要外部因素。

1 分析方法

1.1 Takens定理和相空間重構(gòu)理論

Takens[16]證明了可以找到一個(gè)合適的嵌入維，在這個(gè)嵌入維空間里可以把有規(guī)律的軌跡恢復(fù)出來。該定理表明，可以將單變量時(shí)間序列重構(gòu)成一個(gè)相空間，只要嵌入維數(shù)足夠高，就可構(gòu)建出一個(gè)拓?fù)涞葍r(jià)的相空間。由于混沌系統(tǒng)的策動(dòng)因素是相互影響的，因而在時(shí)間上先后產(chǎn)生的數(shù)據(jù)點(diǎn)也是相關(guān)的，Packard等[17]提出了采用原系統(tǒng)中某變量的延遲坐標(biāo)來重構(gòu)相空間。因此重構(gòu)相空間時(shí)，只對(duì)一個(gè)分量在固定時(shí)延點(diǎn)上進(jìn)行維向量拓展至三維甚至更高維的空間，就可以重構(gòu)一個(gè)拓?fù)涞葍r(jià)的相空間。

設(shè)x(t)，t=1,2,…,N為觀測(cè)時(shí)間序列。重構(gòu)相空間選擇嵌入維數(shù)m，時(shí)間延遲設(shè)為τ，則由x(t)構(gòu)造出一組新的向量序列X(ti)={x(ti),x(ti+τ),…,x(ti+(m-1)τ)}T，i=1,2,…,M，M=N-(m-1)τ，如式(1)所示:

(1)

即

X(t)={X(t1),X(t2),…,X(ti),…,X(tM)}

式(1)表明，時(shí)間延遲τ和嵌入維數(shù)m的選擇是相空間重構(gòu)的主要內(nèi)容，直接影響到相空間重構(gòu)質(zhì)量和其奇異吸引子特征量的描述。

1.2 時(shí)間延遲的選取

對(duì)于無限長(zhǎng)、無噪聲的時(shí)間序列，時(shí)間延遲τ的選擇原則無限制，但對(duì)于有限時(shí)間序列，相空間的特征量則依賴于τ。Taken[16]認(rèn)為時(shí)間延遲和嵌入維數(shù)兩者理論上相互獨(dú)立無關(guān)，兩者可以單獨(dú)確定。也有研究者認(rèn)為兩者的選擇應(yīng)互相關(guān)聯(lián)[18-19]，延遲時(shí)間的選取應(yīng)依賴于延遲時(shí)間窗口，即τw=(m-1)τ。

本文主要采用C-C方法進(jìn)行時(shí)間延遲估計(jì)，并采用互信息算法計(jì)算結(jié)果作為對(duì)比?；バ畔⒘糠椒ㄊ枪烙?jì)重構(gòu)相空間時(shí)間延遲的一種有效方法。該方法利用互信息函數(shù)的第一次極小值來確定相空間重構(gòu)的時(shí)間延遲，F(xiàn)raser等給出了互信息計(jì)算的遞歸算法[20]。Kim等人[19]提出了C-C方法，即通過嵌入時(shí)間序列的關(guān)聯(lián)積分構(gòu)造統(tǒng)計(jì)量來代表非線性時(shí)間序列的相關(guān)性，對(duì)時(shí)間延遲τ和時(shí)間窗口τw同時(shí)進(jìn)行估計(jì)。該方法效果和互信息量法一致，但能有效減少互信息法的計(jì)算量，并能保持其非線性特征，且對(duì)小數(shù)據(jù)組可靠，抗干擾能力較強(qiáng)。

1.3 嵌入維數(shù)的確定

相空間嵌入維數(shù)需要足夠大，以便能刻畫出該系統(tǒng)的奇異吸引子。奇異吸引子是描述混沌系統(tǒng)時(shí)間序列的重要特征，其維數(shù)一般低于相空間的維數(shù)。關(guān)聯(lián)維數(shù)D常用于描述奇異吸引子的維數(shù)，嵌入維數(shù)的選取通常要求m≥2D+1。本文采用飽和關(guān)聯(lián)維法(G-P算法) 求取關(guān)聯(lián)維數(shù)。

Grassberger和Proacaccia[21]依據(jù)是嵌入定理和相空間重構(gòu)理論提出的求取關(guān)聯(lián)維的G-P算法，只需針對(duì)實(shí)驗(yàn)觀測(cè)數(shù)據(jù)就可以得到吸引子維數(shù)，目前采用較多：

對(duì)于嵌入維數(shù)為m的一組向量序列維{Xi}i=1,2,…,N，其中Xi={xi1,xi2,…,xim}為參考點(diǎn)，分別計(jì)算其它N-1個(gè)點(diǎn)與Xi的距離，統(tǒng)計(jì)落在以Xi點(diǎn)為中心且距離小于正數(shù)r的點(diǎn)的個(gè)數(shù)。向量之間的距離一般采用兩種方法，2-范數(shù)和∞-范數(shù)。兩種方法的計(jì)算距離的公式是拓?fù)涞葍r(jià)的，后者計(jì)算量較小[22]。式(2)所示采用的是∞-范數(shù)計(jì)算，即以兩個(gè)向量的最大分量差作為距離:

(2)

滿足該條件的矢量對(duì)數(shù)在種配對(duì)中所占比例被稱為關(guān)聯(lián)積分:

(3)

其中表示Heaviside階躍函數(shù):

(4)

令dm為m維空間中的奇異吸引子最大伸展距離，r≥dm時(shí)Cm(r)=N(N-1)/N2=(N-1)/N,當(dāng)N→∞時(shí)，Cm(r)≈1。實(shí)際上，Cm(r)反映了奇異吸引子中各點(diǎn)之間距離的分布概率，則有:

(5)

根據(jù)分形理論的標(biāo)度不變性，當(dāng)r位于無標(biāo)度區(qū)間時(shí)，D(m,r)是與m和r有關(guān)的常數(shù)，稱為關(guān)聯(lián)指數(shù)，是關(guān)聯(lián)維數(shù)D的確定逼近。取小距離r1和r2，則有：

(6)

對(duì)式6兩邊取對(duì)數(shù)，則有:

(7)

當(dāng)值∣r2-r1∣很小時(shí)，可認(rèn)為D(m,r2) ≈D(m,r1)。則有:

(8)

由式(8)可知，D(m,r)是lnCm(r)—lnr曲線的斜率。當(dāng)r→0時(shí)，可獲得關(guān)聯(lián)維數(shù)D的近似值:

(9)

確定關(guān)聯(lián)維的關(guān)鍵是先確定無標(biāo)度區(qū)間。實(shí)際動(dòng)力系統(tǒng)的分形不同于數(shù)學(xué)上的分形所具有的在無窮尺度上的自相似或自仿射性，而是近似地或統(tǒng)計(jì)意義上存在自相似，這種自相似僅存在于一定的尺度變化范圍，一旦超出了這個(gè)尺度變化范圍，其自相似性就不復(fù)存在，這個(gè)尺度變化范圍就是分形的無標(biāo)度區(qū)。根據(jù)標(biāo)度區(qū)間不變性定理，無標(biāo)度區(qū)間內(nèi)，Cm(r)與r成指數(shù)關(guān)系，如式(9)。標(biāo)量的取值應(yīng)包含在分形的無標(biāo)度區(qū)之內(nèi)，其序列長(zhǎng)度不應(yīng)少于20，以保證估計(jì)到的關(guān)聯(lián)維可靠且滿足相關(guān)性檢驗(yàn)[23]。近些年，針對(duì)無標(biāo)度區(qū)間的確定方法已展開多方面的研究[24- 27]。這些確定的方法有效性和針對(duì)性雖然不相同，但其基本步驟均是：1) 改變嵌入維m值，獲得lnCm(r)-lnr的曲線簇；2) 用最小二乘法擬合每條曲線最佳直線段(無標(biāo)度區(qū))，獲得該直線段的斜率作為嵌入維數(shù)m對(duì)應(yīng)的關(guān)聯(lián)指數(shù)；3) 隨著m的增加，關(guān)聯(lián)指數(shù)會(huì)出現(xiàn)飽和現(xiàn)象，此時(shí)最佳直線段對(duì)應(yīng)的斜率即為該混沌序列對(duì)應(yīng)的關(guān)聯(lián)維數(shù)D。

1.4 最大Lyapunov指數(shù)的估計(jì)方法

混沌系統(tǒng)可通過其吸引子的宏觀特征量表征。宏觀特征量主要包括：關(guān)聯(lián)維，Kolmogorov熵，Lyapunov指數(shù)等。其中，Lyapunov指數(shù)的時(shí)間數(shù)列是根據(jù)序列本身所計(jì)算出來的客觀規(guī)律，是量化初始閉軌道的指數(shù)發(fā)散和估計(jì)方法的混沌量，反映了動(dòng)力系統(tǒng)整體混沌量水平。因此，Lyapunov指數(shù)計(jì)算對(duì)預(yù)測(cè)混沌時(shí)間序列而言顯得尤為重要。

在判斷實(shí)際動(dòng)力系統(tǒng)是否為混沌系統(tǒng)時(shí)，通常只估計(jì)最大Lyapunov指數(shù)λ1。小數(shù)據(jù)量法是只計(jì)算最大Lyapunov指數(shù)的一種方法，是通過對(duì)基本軌道上每個(gè)點(diǎn)的最近鄰近點(diǎn)的平均發(fā)散速率估計(jì)出最大Lyapunov指數(shù)。該算法需要給定嵌入維數(shù)，時(shí)間延遲以及序列的平均周期?？紤]到小數(shù)據(jù)量算法對(duì)嵌入維數(shù)，時(shí)間延遲，數(shù)據(jù)序列長(zhǎng)度以及信噪比的敏感程度較低，具有較好的魯棒性[28]，本文采用了該算法實(shí)現(xiàn)對(duì)最大Lyapunov指數(shù)的求取。但為避免計(jì)算過程中參數(shù)選取的主觀性，采用了C-C方法對(duì)嵌入維數(shù)和時(shí)間延遲同時(shí)估計(jì)。其平均周期的計(jì)算是將序列FFT變換獲得的每個(gè)功率值先與其對(duì)應(yīng)的周期加權(quán)，再最終求加權(quán)平均計(jì)算獲得[29-30]。本文只對(duì)時(shí)間序列的最大Lyapunov指數(shù)進(jìn)行求取，該指數(shù)是否大于零可用來判別一個(gè)時(shí)間序列是否為混沌系統(tǒng)。

1.5 Lyapunov時(shí)間計(jì)算

最大Lyapunov指數(shù)λ1表示的是當(dāng)前兩相鄰軌道距離最大分量的指數(shù)發(fā)散程度：

(10)

式中，δx(0)表示兩相鄰軌道初始距離, δx(t)表示經(jīng)過時(shí)間后的軌道距離。混沌運(yùn)動(dòng)并不是真正的隨機(jī)，他服從確定性的規(guī)律，這表明在一定的臨界距離變化范圍內(nèi)是可以預(yù)測(cè)的，而超過這個(gè)范圍時(shí)，軌道運(yùn)動(dòng)就不可預(yù)測(cè)了[27, 31]。距離變化如式(11)所示，達(dá)到臨界值所經(jīng)歷的時(shí)間如式(12)所示:

(11)

(12)

通常認(rèn)為軌跡分離達(dá)到原間距的數(shù)倍或十幾倍時(shí)，軌道就不確定了。此時(shí)可認(rèn)為lnC≈1，則:

(13)

式(13)所示的時(shí)間為L(zhǎng)yapunov時(shí)間，也稱為最大預(yù)測(cè)時(shí)間。當(dāng)λ1越大，預(yù)測(cè)時(shí)間越短，系統(tǒng)運(yùn)動(dòng)的可預(yù)測(cè)性越差。

2 分析結(jié)果

本文主要對(duì)易北河1997年—2001年的葉綠素a觀測(cè)序列進(jìn)行了混沌特性分析。易北河全長(zhǎng)1094km，為歐洲第四大內(nèi)陸河，流域面積達(dá)148268 km2。流域面積主要覆蓋捷克和德國，葉綠素a的測(cè)點(diǎn)位置為Geesthacht Weir(N53°25.486′ E10°20.665′)，并采用在線葉綠素監(jiān)測(cè)儀(10-AU-005)每小時(shí)采樣1次。本文采用C-C算法[19]進(jìn)行時(shí)間延遲τ的確定，采用G-P算法[21]對(duì)數(shù)據(jù)序列的關(guān)聯(lián)維數(shù)進(jìn)行估計(jì)D，進(jìn)而確定嵌入維數(shù)m，并采用Rosenstein方法[27]計(jì)算了葉綠素a觀測(cè)序列的最大Lyapunov指數(shù)λ1，對(duì)最大預(yù)測(cè)時(shí)間t0進(jìn)行了估計(jì)。

2.1 時(shí)間延遲估計(jì)

每年的葉綠素a觀測(cè)數(shù)據(jù)序列通過C-C方法獲得的時(shí)間延遲(圖1)。時(shí)間延遲τ=tτs，τs為采樣時(shí)間間隔，取值為1，則τ=t。曲線第一個(gè)極小值所對(duì)應(yīng)的t為時(shí)間延遲，圖中所示時(shí)間延遲為8或12，該遲延與互信息法所確定的時(shí)間基本一致(圖1)。互信息法中時(shí)間延遲對(duì)應(yīng)曲線第一個(gè)極小值。

時(shí)間延遲τ比較如表1所示。采用互信息方法和C-C方法獲得τ差異較小，表明這兩種方法在τ的判斷上比較接近。本文選用C-C方法獲得的τ值進(jìn)行關(guān)聯(lián)維數(shù)的計(jì)算。

2.2 關(guān)聯(lián)維數(shù)估計(jì)

針對(duì)1997—2001年的葉綠素a觀測(cè)序列采取了C-C方法獲得時(shí)間延遲τ (表1)，并采用G-P算法獲取了對(duì)應(yīng)各年份的關(guān)聯(lián)維數(shù)D，如圖2所示，其中標(biāo)度r起始值分別參考空間向量距離的最小值和最大值，其步長(zhǎng)設(shè)定為30。無標(biāo)度區(qū)間的確定方法是采用目前常用的視覺識(shí)別方法。

圖2所示，每年的葉綠素a的觀測(cè)序列均具有低維分形維。通過C-C算法獲得的最佳嵌入維數(shù)與通過G-P算法中獲得的飽和嵌入維數(shù)相接近，如表2所示。飽和嵌入維數(shù)是指最佳直線段斜率剛達(dá)到飽和時(shí)所對(duì)應(yīng)的嵌入維數(shù)。當(dāng)增加超過飽和嵌入維數(shù)時(shí)，其直線段斜率變化趨于穩(wěn)定，由G-P算法可知，對(duì)應(yīng)的最佳直線段斜率值即為關(guān)聯(lián)維數(shù)。

圖1 葉綠素a觀測(cè)序列時(shí)間延遲τ (τ=t)Fig.1 The time delay of the sequence of hourly chlorophyll a observation (τ=t)

表1 葉綠素a觀測(cè)序列時(shí)間延遲不同算法計(jì)算結(jié)果比較*

表2 各年關(guān)聯(lián)維數(shù)和嵌入維數(shù)列表

2.3 最大Lyapunov指數(shù)估計(jì)

選取C-C方法計(jì)算所得的嵌入維數(shù)m和時(shí)間延遲τ，如表1，根據(jù)Rosenstein算法[28]，結(jié)果如圖3所示。圖中橫坐標(biāo)表示為離散時(shí)間步長(zhǎng)，縱坐標(biāo)ln表征近鄰點(diǎn)對(duì)在第i個(gè)離散時(shí)間步長(zhǎng)后的距離對(duì)數(shù)值。圖3所示，在i=20—100的范圍內(nèi)，曲線表現(xiàn)出較好的直線特性。選取該區(qū)域，采用最小二乘法做出回歸直線，該直線的斜率即為最大Lyapunov指數(shù)λ1。Rosenstein算法中根據(jù)時(shí)間序列的平均周期限制相空間中最近鄰點(diǎn)的短暫分離，時(shí)間序列平均周期的算法參考Rathje等[29-30]提出的算法。每年葉綠素a濃度時(shí)間序列的平均周期如圖3中tp值所示。

如表3所示，各年的葉綠素a觀測(cè)序列的最大Lyapunov指數(shù)均為非零數(shù)，表明其時(shí)間序列具有混沌特性，Lyapunov指數(shù)接近于零值，表明該時(shí)間序列混沌特性較弱。5a的葉綠素a觀測(cè)序列最大預(yù)測(cè)時(shí)間最大約為19 d，最短約為8 d，平均約為14 d (兩周)。

表3 不同年份最大Lyapunov指數(shù)和最大預(yù)測(cè)時(shí)間列表

圖2 葉綠素a觀測(cè)序列關(guān)聯(lián)維數(shù)Fig.2 The correlation dimension D of hourly chlorophyll a observation sequence

圖3 不同年份葉綠素a觀測(cè)序列Lyapunov指數(shù)(Rosenstein算法)Fig.3 The largest Lyapunov index of the chlorophyll a sequence estimated by the small-data method

2.4 河流日流量時(shí)間序列混沌特性分析

采用與葉綠素a觀測(cè)序列同樣的分析方法，對(duì)易北河1997年3月3日至2001年10月31日之間的日徑流量時(shí)間序列進(jìn)行混沌特性分析，流量采樣點(diǎn)位置為Neu Darcha((N53°14.207′ E10°53.314′)，采樣設(shè)備為多普勒測(cè)流儀，每日采樣一次。首先采用了C-C方法同時(shí)對(duì)延遲時(shí)間窗口和時(shí)間延遲同時(shí)估計(jì)。C-C方法中的最小值對(duì)應(yīng)時(shí)間延遲窗口值，的第一個(gè)零點(diǎn)或者的第一個(gè)極小值對(duì)應(yīng)著時(shí)間延遲(τ=t)。由圖4可知，tw=60，τ=17，其中的第一個(gè)零點(diǎn)對(duì)應(yīng)值為τ=60。通常取較小值作為時(shí)間延遲，避免使得非線性模型的擬合關(guān)系復(fù)雜，降低難度和減少計(jì)算量[24]。C-C方法后的時(shí)間延遲與互信息法計(jì)算結(jié)果相接近，后者τ=15，對(duì)應(yīng)曲線的第一個(gè)極小值。

圖4 1997—2001年易北河日徑流量時(shí)間序列時(shí)間延遲Fig.4 The time delay of the daily runoff sequence (1997—2001)

1997—2001年易北河日徑流量時(shí)間序列飽和關(guān)聯(lián)維數(shù)為9，如圖5所示。而C-C方法獲得的tw和τ，根據(jù)可得tw=(m-1)τ，m=5。G-P算法獲得的關(guān)聯(lián)維數(shù)D=1.84，則日徑流量時(shí)間序列的嵌入維數(shù)m≥2D+1。同樣，C-C算法獲得的嵌入維數(shù)也滿足此條件。

圖5 1997—2001年易北河日徑流量時(shí)間序列關(guān)聯(lián)維數(shù)Fig.5 The correlation dimension of the daily runoff sequence (1997—2001)

根據(jù)Rosenstein算法，取m=5，τ=17，tp=26(平均周期)，計(jì)算結(jié)果如圖6，由圖中觀察可知，在0—100的范圍內(nèi)具有直線上升趨勢(shì)，采用最小二乘法對(duì)該范圍進(jìn)行直線擬合，獲得直線段斜率，即估計(jì)最大Lyapunov指數(shù)為λ1=0.0125，估計(jì)最大預(yù)測(cè)時(shí)間為80d。非零數(shù)的指數(shù)估計(jì)和關(guān)聯(lián)維數(shù)表明日徑流量同樣具有混沌特性。

圖6 1997—2001年易北河日徑流量時(shí)間序列Lyapunov指數(shù)(Rosenstein算法)Fig.6 The largest Lyapunov index of daily runoff sequence (1997—2001)

3 結(jié)果討論

采用兩種方法即C-C方法和互信息法進(jìn)行了時(shí)間序列時(shí)間延遲的判定，從葉綠素a觀測(cè)序列和日徑流量時(shí)間序列的分析結(jié)果可知，兩種方法各自獲得的時(shí)間延遲比較接近，如表1，圖4所示。在重構(gòu)相空間時(shí)，τ不取最佳值只會(huì)影響關(guān)聯(lián)維數(shù)的計(jì)算，并不會(huì)影響重構(gòu)吸引子如實(shí)反映系統(tǒng)的動(dòng)力學(xué)性質(zhì)[27]。τ取值太小，迫使空間軌跡沿同一方向擠壓，信息被壓縮，獲取困難。反之，導(dǎo)致相鄰時(shí)刻的動(dòng)力學(xué)狀態(tài)變化劇烈，使得簡(jiǎn)單的系統(tǒng)運(yùn)動(dòng)變得看起來很復(fù)雜。采用C-C方法獲得的τ依賴于時(shí)間延遲窗口τw，τw反映了數(shù)據(jù)依賴的最大時(shí)間，而互信息方法僅僅反映了第一次局部最大時(shí)間，不能估計(jì)τw。另外，評(píng)價(jià)重構(gòu)吸引子的質(zhì)量是其幾何上的重復(fù)性和不相關(guān)性，重構(gòu)吸引子應(yīng)具有較低的重復(fù)性和較強(qiáng)的相關(guān)性。因此，τw是一種用于估計(jì)維數(shù)的更好的量。盡管C-C方法是利用統(tǒng)計(jì)結(jié)果獲得，但通過和互信息對(duì)比表明，在確定序列時(shí)間延遲上仍具有較好的效果，而且該方法能同時(shí)確定嵌入窗口和嵌入維數(shù)，為較好的估計(jì)Lyapunov指數(shù)提供了前提條件。

為更好的驗(yàn)證C-C方法計(jì)算嵌入維數(shù)的可靠性，采用了G-P算法對(duì)嵌入維數(shù)進(jìn)行估計(jì)。計(jì)算結(jié)果表明，通過G-P算法獲得的飽和嵌入維數(shù)與C-C方法獲得嵌入維數(shù)相近，如圖2和圖5所示。這進(jìn)一步驗(yàn)證了C-C方法的可靠性。

根據(jù)G-P算法獲得關(guān)聯(lián)維數(shù)過程中，對(duì)無標(biāo)度區(qū)的估計(jì)采用了視覺識(shí)別法。無標(biāo)度區(qū)的判斷盡管有很多種方法[25-26, 32]，但在無法獲知理論關(guān)聯(lián)維數(shù)的情況下，視覺識(shí)別法仍可作為一種方便快捷的估計(jì)方法。如圖2所示，當(dāng)m取較小值時(shí)，其直線段(無標(biāo)度區(qū))較為明顯，直線區(qū)間范圍較大。隨著m的增加，其直線段范圍發(fā)生變化，并有逐漸變小的趨勢(shì)，如圖2中的曲線包圍范圍所示，這與文獻(xiàn)[23]描述結(jié)果一致。該方法獲得關(guān)聯(lián)維數(shù)的精度對(duì)判定系統(tǒng)混沌特性以及確定關(guān)聯(lián)維數(shù)并無顯著影響[33]，而且由C-C方法確定的嵌入維數(shù)均比獲得由關(guān)聯(lián)維數(shù)確定的最小嵌入維數(shù)大得多(>2D+1)，確保展開的相空間與原系統(tǒng)的動(dòng)力學(xué)特性拓?fù)涞葍r(jià)。同樣，葉綠素a序列的低關(guān)聯(lián)維數(shù)估計(jì)值表明藻類生長(zhǎng)在較低維空間產(chǎn)生自相似，但藻類生長(zhǎng)的復(fù)雜性可能使得需要采取高的嵌入維數(shù)值才能完全展現(xiàn)藻類的生長(zhǎng)軌跡特征。

采用Rosenstein算法對(duì)最大Lyapunov指數(shù)進(jìn)行了估計(jì)，先后對(duì)葉綠素a和日徑流量時(shí)間序列進(jìn)行了指數(shù)估計(jì)。結(jié)果表明，兩種時(shí)間序列的最大Lyapunov指數(shù)估計(jì)值均較小，如表2，圖6所示，這表明該序列的鄰近吸引子軌道離散程度較低，可以進(jìn)行較長(zhǎng)時(shí)期的預(yù)測(cè)。其中，5a的日凈流量的最大預(yù)測(cè)時(shí)間估計(jì)為80d，5a的葉綠素a最大預(yù)測(cè)時(shí)間平均為14d，前者要遠(yuǎn)大于后者?？梢?，盡管河流徑流量對(duì)葉綠素a有影響，但可能并不是對(duì)后者產(chǎn)生混沌特性影響的主要因素。

河流水體中的藻中生長(zhǎng)的另一個(gè)重要因素是天氣條件，而氣候系統(tǒng)已被證實(shí)是具有混沌特性。目前天氣預(yù)報(bào)可以被較好預(yù)測(cè)的時(shí)間約為15d，和葉綠素a觀測(cè)序列的最大預(yù)測(cè)時(shí)間相吻合。相比于水文條件，氣候條件可能是對(duì)葉綠素a的混沌特性產(chǎn)生影響的一個(gè)主要因素。

在Elbe河流中的硅藻幾乎總是優(yōu)勢(shì)藻種，其限制性營養(yǎng)鹽為硅。因此，葉綠素a所反映的主要是硅藻濃度，始終作為優(yōu)勢(shì)藻，使得由于藻種間的競(jìng)爭(zhēng)導(dǎo)致的混沌現(xiàn)象似乎不可能發(fā)生，另外，其限制性營養(yǎng)鹽硅輸入主要是來自于Elbe河源頭的Elbe砂巖山的硅砂巖土層，而沿途補(bǔ)充的營養(yǎng)源多為氮磷，因此，由于營養(yǎng)鹽輸入變化導(dǎo)致的藻種濃度出現(xiàn)混沌特征可以忽略。這也說明了氣象因素要相對(duì)于其他因素更有可能造成藻種葉綠素a濃度的混沌變化。

4 總結(jié)

采用了混沌理論對(duì)易北河葉綠素a1997年至2001年間各年的觀測(cè)序列進(jìn)行了混沌分析，分析結(jié)果表明，各年的葉綠素a觀測(cè)序列均具有低維混沌特性，同樣在對(duì)5a易北河連續(xù)日徑流量時(shí)間序列分析表明，該徑流量也具有低維混沌特性。通過對(duì)兩種序列的最大可預(yù)測(cè)時(shí)間比對(duì)表明徑流量等水文因子可能并不是影響葉綠素a混沌特性的主要因素，相反氣象因子可能占據(jù)主導(dǎo)地位。

致謝：Mr.Scharfe博士提供德國易北河歷史數(shù)據(jù),特此致謝。

[1] Huisman J, Weissing F J. Biodiversity of plankton by species oscillations and chaos. Nature, 1999, 402(6760): 407- 410.

[2] Connell J H. Diversity in tropical rain forests and coral reefs. Science, 1978, 199(4335): 1302- 1310.

[3] Sala S, Vighi M. GIS-based procedure for site-specific risk assessment of pesticides for aquatic ecosystems. Ecotoxicology and Environmental Safety, 2008, 69(1): 1- 12.

[4] Malmaeus J M, Blenckner T, Markensten H, Persson I. Lake phosphorus dynamics and climate warming: a mechanistic model approach. Ecological Modelling, 2006, 190(1/2): 1- 14.

[5] Hu W P, J?rgensen S E, Zhang F B. A vertical-compressed three-dimensional ecological model in Lake Taihu, China. Ecological Modelling, 2006, 190(3/4): 367- 398.

[6] 饒群, 芮孝芳. 完全混合系統(tǒng)總磷隨機(jī)模型研究. 水科學(xué)進(jìn)展, 2002, 13(1): 21- 25.

[7] Xu F L, J?rgensen S E, Tao S, Li B G. Modeling the effects of ecological engineering on ecosystem health of a shallow eutrophic Chinese lake (Lake Chao). Ecological Modelling, 1999, 117(2/3): 239- 260.

[8] Hassan H, Hanaki K, Matsuo T. A modeling approach to simulate impact of climate change in lake water quality: phytoplankton growth rate assessment. Water Science and Technology, 1998, 37(2): 177- 185.

[9] 趙曉東, 潘江, 李金頁, 陶曉磊, 龐坤. 銅綠微囊藻和斜生柵藻非穩(wěn)態(tài)營養(yǎng)鹽限制條件下的生長(zhǎng)競(jìng)爭(zhēng)特性. 生態(tài)學(xué)報(bào), 2011, 31(13): 3710- 3719.

[10] Scharfe M, Callies U, Bl?cker G, Petersen W, Schroeder F. A simple Lagrangian model to simulate temporal variability of algae in the Elbe River. Ecological Modelling, 2009, 220(18): 2173- 2186.

[11] Tsonis A A, Elsner J B. The weather attractor over very short timescales. Nature, 1988, 333(6173): 545- 547.

[12] Radhakrishnan P, Dinesh S. An alternative approach to characterize time series data: case study on Malaysian rainfall data. Chaos, Solitons and Fractals, 2006, 27(2): 511- 518.

[13] 丁晶, 王文圣, 趙永龍. 長(zhǎng)江日流量混沌變化特性研究——Ⅱ相空間嵌入維數(shù)的確定. 水科學(xué)進(jìn)展, 2003, 14(4): 412- 416.

[14] 丁晶, 王文圣, 趙永龍. 長(zhǎng)江日流量混沌變化特性研究——Ⅰ相空間嵌入滯時(shí)的確定. 水科學(xué)進(jìn)展, 2003, 14(4): 407- 411.

[15] Sivakumar B. Rainfall dynamics at different temporal scales: a chaotic perspective. Hydrology and Earth System Sciences, 2001, 5(4): 645- 652.

[16] Takens F. Detecting strange attractors in turbulence. Lecture Notes in Mathematics, 1981, 898: 366- 381.

[17] Packard N H, Crutchfield J P, Farmer J D, Shaw R S. Geometry from a time series. Physical Review Letters, 1980, 45(9): 712- 716.

[18] Kugiumtzis D. State space reconstruction parameters in the analysis of chaotic time series-the role of the time window length. Physica D: Nonlinear Phenomena, 1996, 95(1): 13- 28.

[19] Kim H S, Eykholt R, Salas J D. Nonlinear dynamics, delay times, and embedding windows. Physica D: Nonlinear Phenomena, 1999, 127(1/2): 48- 60.

[20] Fraser A M, Swinney H L. Independent coordinates for strange attractors from mutual information. Physical Review A, 1986, 33(2): 1134- 1140.

[21] Grassberger P, Procaccia I. Characterization of strange attractors. Physical Review Letters, 1983, 50(5): 346- 349.

[22] 王安良, 楊春信. 評(píng)價(jià)奇怪吸引子分形特征的Grassberger-Procaccia算法. 物理學(xué)報(bào), 2002, 51(12): 2719- 2729.

[23] 黨建武, 黃建國. 基于G. P算法的關(guān)聯(lián)維計(jì)算中參數(shù)取值的研究. 計(jì)算機(jī)應(yīng)用研究, 2004, 21(1): 48- 51.

[24] 韓敏. 混沌時(shí)間序列預(yù)測(cè)理論與方法. 北京: 中國水利水電出版社, 2007.

[25] 黨建武, 施怡, 黃建國. 分形研究中無標(biāo)度區(qū)的計(jì)算機(jī)識(shí)別. 計(jì)算機(jī)工程與應(yīng)用, 2003, 39(12): 25- 27.

[26] 巫兆聰. 分形分析中的無標(biāo)度區(qū)確定問題. 測(cè)繪學(xué)報(bào), 2002, 31(3): 240- 244.

[27] 呂金虎, 陸君安, 陳士華. 混沌時(shí)間序列分析及其應(yīng)用. 武漢: 武漢大學(xué)出版社, 2002.

[28] Rosenstein M T, Collins J J, De Luca C J. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena, 1993, 65(1/2): 117- 134.

[29] Rathje E M, Abrahamson N A, Bray J D. Simplified frequency content estimates of earthquake ground motions. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(2): 150- 159.

[30] Rathje E M, Faraj F, Russell S, Bray J D. Empirical relationships for frequency content parameters of earthquake ground motions. Earthquake Spectra, 2004, 20(1): 119- 144.

[31] 劉式達(dá), 梁福明, 劉式適, 辛國君. 自然科學(xué)中的混沌和分形. 北京: 北京大學(xué)出版社, 2003.

[32] 陳敏. 一種簡(jiǎn)單的GP算法無標(biāo)度區(qū)識(shí)別方法. 計(jì)算機(jī)與信息技術(shù), 2007, (11): 35- 36, 39- 39.

[33] 姬翠翠, 朱華, 江煒. 混沌時(shí)間序列關(guān)聯(lián)維數(shù)計(jì)算中無標(biāo)度區(qū)間識(shí)別的新方法. 科學(xué)通報(bào), 2010, 55(31): 3069- 3076.

[34] Smith L A. Intrinsic limits on dimension calculations. Physics Letters A, 1988, 133(6): 283- 288.

[35] 洪時(shí)明, 洪時(shí)中. 用Grassberger-Procaccia方法計(jì)算吸引子維數(shù)的基本限制. 物理學(xué)報(bào), 1994, 43(8): 1228- 1233.

[36] 林振山. 長(zhǎng)期預(yù)報(bào)的相空間理論和模式. 北京: 氣象出版社, 1993.

Analysis of the growth characteristics of river algae based on chaos theory: a case study of Elbe River

ZHAO Xiaodong1,2,*, ZHANG Hongjian1, ZHOU Hongliang1

1StateKeyLaboratoryofIndustrialControlTechnology,ZhejiangUniversity,Hangzhou310007,China2CollegeofMetrologicalandMeasurementEngineering,ChinaJiliangUniversity,Hangzhou310018,China

Algal growth at equilibrium is not sustainable in a river or lake experiencing point or non-point source pollution, which is likely to change with season, location, and human activity. Chlorophyll a, as a common indictor, is an important reference point for water resource management. It can be affected by the biodiversity of species, as a result of their oscillations, and chaotic fluctuations. Algal growth has many highly nonlinear characteristics. The characteristics that differ among species, however, are all susceptible to change under external disturbance. It is difficult to provide a comprehensive and detailed description of nonlinear algal growth. Variations in chlorophyll a tend to maintain certain regularity; for example, seasonal variation and the 24 hour cycle, which also display self-similarity. However, it is difficult to observe similar variations at different times from the studied sampling series. These features are similar to aspects of chaotic motion, such as boundedness, ergodicity, and inherent randomness. The sampling of the indictor chlorophyll a is typically performed on an hourly, weekly, or even monthly basis. With higher sampling frequency, the sampling series of chlorophyllain the field becomes more unstable and appears to be more chaotic. Therefore, this paper aims to study the variations in a chlorophyll a series sampled from the field, rather than constructing a theoretical model to recapitulate the field data. While there is extensive research on algae in many lakes or rivers, few studies discuss the prediction of algal growth times in aquatic environments. The aspects of algal growth times are partially addressed in this study.The variation characteristics of the algal data series were analyzed using chaos theory. The characters of the one-dimensional time series were recovered by reconstruction into the multi-dimensional phase space, using phase space reconstruction. The reconstruction parameters, namely the embedding dimensionmand time delayτ, were estimated using the correlation integral method (C-C method). The correlation dimension,D, is the basic mathematical description of the strange attractor, which is the main characteristic of a chaotic system.Dwas calculated using the Grassberger-Procaccia algorithm (G-P algorithm). Only the largest Lyapunov exponent,λ1, was estimated through the small data method to evaluate the diffusion degree of the phase trajectory. The reciprocal ofλ1is the upper bound of the deterministic prediction time in the chaotic system, which is designated as the Lyapunov timet0. This property indicates that the system is unpredictable beyondt0. In this paper, the chaotic characteristics of hourly chlorophyll a concentrations and the daily runoff time series of the Elbe River over a five year period (1997—2001) were analyzed. It was found that both sequences had the properties of low-dimension chaos withλ1>0,D= 2.75—4.02 for the chlorophyll a series, andD= 1.84 for the runoff series. The average value oft0(14 days) was estimated for the five-year chlorophyll a sequence data. These findings were remarkably close to the current biggest day-to-day weather forecast time (two to three weeks). A much longer value oft0for the runoff series, for the same period, was 80 days. This result indicated that compared to the weather factors, the runoff factor is clearly weaker in affecting the chaotic characteristics of chlorophyll a.

algal growth; chaos; correlation dimension; Lyapunov exponent

中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助(2014FZA5008)

2013- 11- 05;

日期:2014- 11- 03

10.5846/stxb201311052669

*通訊作者Corresponding author.E-mail: zhaoxd@cjlu.edu.cn

趙曉東, 張宏建, 周洪亮.基于混沌理論的河流藻類生長(zhǎng)特性分析——以德國易北河為例.生態(tài)學(xué)報(bào),2015,35(17):5585- 5596.

Zhao X D, Zhang H J, Zhou H L.Analysis of the growth characteristics of river algae based on chaos theory: a case study of Elbe River.Acta Ecologica Sinica,2015,35(17):5585- 5596.