趙曉東, 張宏建, 周洪亮

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

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

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

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

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

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

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

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

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

1 分析方法

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

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

設(shè)x(t)，t=1,2,…,N為觀測時間序列。重構(gòu)相空間選擇嵌入維數(shù)m，時間延遲設(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ù)m的選擇是相空間重構(gòu)的主要內(nèi)容，直接影響到相空間重構(gòu)質(zhì)量和其奇異吸引子特征量的描述。

1.2 時間延遲的選取

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

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

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

相空間嵌入維數(shù)需要足夠大，以便能刻畫出該系統(tǒng)的奇異吸引子。奇異吸引子是描述混沌系統(tǒng)時間序列的重要特征，其維數(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算法，只需針對實驗觀測數(shù)據(jù)就可以得到吸引子維數(shù)，目前采用較多：

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

(2)

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

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

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

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

1.5 Lyapunov時間計算

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

(10)

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

(11)

(12)

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

(13)

式(13)所示的時間為Lyapunov時間，也稱為最大預(yù)測時間。當(dāng)λ1越大，預(yù)測時間越短，系統(tǒng)運動的可預(yù)測性越差。

2 分析結(jié)果

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

2.1 時間延遲估計

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4 河流日流量時間序列混沌特性分析

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

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

1997—2001年易北河日徑流量時間序列飽和關(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ù)m≥2D+1。同樣，C-C算法獲得的嵌入維數(shù)也滿足此條件。

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

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

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

3 結(jié)果討論

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

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

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

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

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

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

4 總結(jié)

采用了混沌理論對易北河葉綠素a1997年至2001年間各年的觀測序列進行了混沌分析，分析結(jié)果表明，各年的葉綠素a觀測序列均具有低維混沌特性，同樣在對5a易北河連續(xù)日徑流量時間序列分析表明，該徑流量也具有低維混沌特性。通過對兩種序列的最大可預(yù)測時間比對表明徑流量等水文因子可能并不是影響葉綠素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)總磷隨機模型研究. 水科學(xué)進展, 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)鹽限制條件下的生長競爭特性. 生態(tài)學(xué)報, 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] 丁晶, 王文圣, 趙永龍. 長江日流量混沌變化特性研究——Ⅱ相空間嵌入維數(shù)的確定. 水科學(xué)進展, 2003, 14(4): 412- 416.

[14] 丁晶, 王文圣, 趙永龍. 長江日流量混沌變化特性研究——Ⅰ相空間嵌入滯時的確定. 水科學(xué)進展, 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] 王安良, 楊春信. 評價奇怪吸引子分形特征的Grassberger-Procaccia算法. 物理學(xué)報, 2002, 51(12): 2719- 2729.

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

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

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

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

[27] 呂金虎, 陸君安, 陳士華. 混沌時間序列分析及其應(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] 劉式達, 梁福明, 劉式適, 辛國君. 自然科學(xué)中的混沌和分形. 北京: 北京大學(xué)出版社, 2003.

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

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

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

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

[36] 林振山. 長期預(yù)報的相空間理論和模式. 北京: 氣象出版社, 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

中央高校基本科研業(yè)務(wù)費專項資金資助(2014FZA5008)

2013- 11- 05;

日期:2014- 11- 03

10.5846/stxb201311052669

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

趙曉東, 張宏建, 周洪亮.基于混沌理論的河流藻類生長特性分析——以德國易北河為例.生態(tài)學(xué)報,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.