復雜時間序列的多尺度分析1)

高劍波

(北京師范大學地理科學學部地理數據與應用分析中心,北京 100875)

引言

20 世紀科學巨匠馮·諾伊曼臨謝世時,對死亡后的人不再有思維,感到迷惑和恐懼.自2021 年8 月25 日吾師鄭哲敏先生逝世至今,我有時也有這樣的感覺,因不能與先生電話或面談而得到先生對某件事的看法.但有時又覺得好像不全如此,因為先生的思維方式至少部分傳下來了.在這篇紀念先生的文章里,我想在談學術研究前,先提一下先生對我個人的影響.這些影響基本決定了我在非線性領域和其他領域研究問題的偏好.本文主要介紹先生和我在非線性領域的一些工作,及后來這些工作如何被拓展并演變成一些通用的非線性問題的分析方法.這些方法已被廣泛用于自然科學、工程技術和社會科學的諸多領域.它們特別適用于各領域(包括運維)的故障診斷,生物醫學數據的分析,及不確定性的度量.大家要是覺得有些方法對研究有用,請隨時聯系我,索要相關程序和使用說明.

1 教誨

1988 年,我從浙江大學工業自動化專業畢業.在浙大時,跟數學老師周先意先生,學了不少數學.周老師是陳建功和蘇步青先生的學生,是數學家王元的同班同學.周老師在1958 年被發配到江西農村.空余時間都在思考教數學時如何把分析、代數和幾何融合起來.1978 年蘇步青先生將其調回浙大.由周老師教我們數學,是我們的幸運.周老師對學生極好,跟我們聊天時經常提到4 個他最得意的學生,曹銀和博士便是其中之一,其時他已在北京力學所,師從白以龍院士,跟鄭先生也熟.來京前,我沒見過曹師兄,但經常與其通信交流.

1987 年大三時,在浙大電機系,我聽了一個外校教授關于混沌動力學的報告,覺得這個領域太迷人了,因此跟曹師兄提及,并表示以后很想做這方面的研究.其時力學所正籌劃于次年6 月建立非線性連續介質力學開放實驗室(LNM),并由鄭先生掛帥,所以曹師兄問我是否愿意去力學所跟鄭先生讀研.我自然很激動.之后曹師兄又帶來一個好消息,說鄭先生會幫我辦理保送手續,這樣就不用擔心因沒學過力學而考不上(不像現在考研,當時只考專業課),這些溝通發生在1987 年10 月.至11 月下旬,出了些問題,曹師兄跟我說,鄭先生因太忙,忘了去給我申請免試讀研,而申請的截止期是11 月1 號.所以師兄問我,是否等一年再說.我想了想,說還是考吧.考試在2 月初.我花了2 個月,學了4 門力學課,考上了力學所.

我在力學所待了6 年.第一年在玉泉路的研究生院,與先生幾乎沒有接觸.后5 年,幾乎隨時可以找先生聊.在LNM 剛開始時,真正做非線性的年輕人就兩三位,包括何國威院士、趙力師兄和我.可惜,我因資質有限,當時學得很有限.下面我講五個方面.

第一點與情懷有關.在追思先生時,很多人提到先生說過的一句話,“不要做錦上添花的事,要做雪中送炭的事”.這個情懷的精髓是堅持不懈地努力,做出最好的工作.這包含兩個方面.一是對自己的要求.聽說先生90 歲高齡時還自己推導公式.二是幫助學生和年輕人,做出最好的工作.記得2007 年,我寄給先生我的英文著作[1]后,先生很喜歡,隨即建議我寄一本給馮元禎先生.馮先生是生物醫學工程奠基人之一,生物力學之父,是鄭先生最好的朋友.馮先生當時已是88 歲高齡.開會時經常碰到加州大學的教授說,馮先生如此高齡,每年還會在PNAS 上發表高質量的文章,太令人敬佩.在書寄給馮先生后,近一個月,沒收到任何回音,我忐忑不安,心里默念是不是自己書寫得不好或馮先生的專業與我相差較多.此時,馮先生的親筆書信到了,信中特別強調,他在每一頁都學到了一些新東西.我特別激動,告知了鄭先生信的內容,先生也很高興.馮先生的鼓勵對我影響很大.馮先生在信中還建議我看看他的好友,美國工程院院士黃鍔(Norden Huang)關于經驗模態分解(empirical mode decomposition,EMD)的工作[2].馮先生太有洞見了！他看出了我們數據分析的體系還缺類似于EMD 的一塊.所幸,當時我們已經發展了一個新方法.這在第4 節會講到.總的來說,馮先生的提醒,很大地促進了我們發展自適應去趨勢、去噪、多尺度分解和分形分析.

另一件事,好像是2010 年初秋,先生還在美國探訪兒子.我回國探親,回美時,路過北京.遠在美國的先生安排我去拜訪林家翹先生.因郭永懷先生和李佩先生是林先生的好友,鄭先生托李佩先生帶我去清華拜訪林先生.當時林先生已是95 歲的高齡了.到了林先生的別墅,聊了一個多小時,只聊一個主題,如何理解蛋白質結構.林先生最出名的工作是流體力學的流動穩定性和星系螺旋結構的密度波理論,但關于此卻一字未提.不囿于自己熟悉的領域,隨著時代發展不斷拓展新的領域,這就是大師們的情懷.身處百年未遇之大變局時期,我們必須發揚光大大師們的這種求真務實和開拓創新的精神.

第二點與心態有關.2002 年8 月,我在佛羅里達大學當助理教授,其時離開力學所和鄭先生已有8 年,自己首次帶研究生,我想,怎么教能好些呢？腦袋里突然冒出來一句先生跟我講過的話:100%的相信,同時又要120%的懷疑.先生的高明之處在于,不用另外解釋,我自然覺得這句話很有道理.在我指導學生時,還展開解釋了,所以落了下乘:前半句,100%的相信,表示對所選的領域和問題的重要性和前瞻性沒有任何懷疑,因此會全身心、充滿激情地投入到學習和研究中,只有這樣,才能將天賦真正變成才能,并發揮出來.后半句,表示在所選的領域中,一些最基本的東西都可能有錯.說的俗一些,就是不要被權威所束縛.其實,樸實的心態應該是,充分欣賞領域里那些漂亮的工作,但不該有誰是權威、能否超越等想法,而應該努力找到所選問題的最佳解,并成為最好的專家之一.雖然我不如先生高明,但這個教學方法還是有效的,如我在佛羅里達大學的博士生胡靜,深得此法的精髓,對后面第3，4 節描述的方法,有很大的貢獻.

第三點與思考有關.先生與我聊天時,強調最多的是,要自由思考.剛到力學所時,有朋友勸我一定要多讀鄭先生出名的那些文章.但因我缺乏基礎,同時當時看不到那些工作與非線性科學的關系,且鄭先生也沒有要求我讀,所以在力學所的幾年,基本沒讀過鄭先生的文章.2014 年,在慶祝鄭先生90 歲生日的會議上,與我很熟的一位老教授,談慶明老師,跟我說,當年我進力學所時,鄭先生特意囑咐他們,不要限制小高,讓他完全自由發展.談老師特別強調,鄭先生這樣關照的,我是絕無僅有的一個.聽談老師這么說,我很感動,也很慚愧.我想主要是我沒力學背景,所以鄭先生特別關照我.當然,由此也可看出鄭先生對非線性科學的重視.可惜,1994 年出國后,有近10 年時間,我只斷斷續續地作為業余愛好者繼續做些非線性問題的研究,所以成績很有限.

第四點是寬容.事情得從1989 年秋天我在LNM 做的報告說起.鄭先生請了很多人來聽,包括北京大學的朱照宣先生.那是我做的第一個報告,講得自然很不好,很多人聽后有些不滿意.不過,會上鄭先生沒有批評我,會后更沒約束我的思維.

后來我沒有轉行做量子計算,理由是,大尺度、宏觀的量子糾纏態幾十年內是沒法實現的,因此不可能發展出能用于通用計算的量子計算機.這個推理不一定嚴謹,但決定是正確的,21 世紀量子計算、量子通信需要的是實驗方面的俊才,但于此我不擅長.另外,1993 年,MIT 的Cuomo 和 Oppenheim[3-4]用Pecora 和 Carroll[5]發展的混沌同步開創了用混沌做安全通訊的新領域.在1998 年前后,美國海軍研究部的Pecora 博士組織了一個多學科大學研究計劃(Multidisciplinary University Research Initiative,MURI,類似 DARPA 的項目) 以推動混沌安全通訊的軍用和民用的研究.參加項目的有加州大學圣地亞哥分校(UCSD)、加州大學洛杉磯分校(UCLA)和斯坦福大學.其時,我在UCLA 電子工程系跟隨Rubin 教授讀博.系里的兩位老師,Jia-Ming Liu 教授和Kung Yao 教授,參與了MURI 項目.姚教授做通訊,劉教授用半導體激光實現通訊.我對此很感興趣,參與了劉教授團隊的工作,因此與幾位同學,Steven Huang、Howard Chen、Shuo Tang 和 C.C.Chen 成了好朋友.前三位是劉教授的學生,最后一位是姚教授的學生.除了Shuo Tang 去英屬哥倫比亞大學當教授,其他三位都去了臺灣當教授.與他們討論時,我告知C.C.Chen,混沌通訊的功耗應該比Cuomo等[5]所說的大,因Cuomo 等在模擬隨機微分方程時,把隨機項處理錯了(成正比,但Cuomo 等用dt代替了 dt—當 dt?1 時,?dt,相當于將實際噪聲縮小了幾個數量級,因而得出功耗小的結論).C.C.Chen 和姚教授發現我是對的,因而證明,混沌通訊也許能用于軍用,但在民用領域,價值很有限[6-7].之后,這方面研究的熱度在美國大降,但在香港和歐洲仍持續了很久.

第五點一定要重視科學實驗和實驗數據分析,盡量從重大的現實問題中,抽象出具有根本科學意義的問題,然后同時解決現實和科學問題.我想這是從普朗特傳至馮·卡門、馮·卡門傳至錢學森、再由錢學森傳至鄭哲敏的技術科學的精髓.

2 混沌時間序列分析

為了不讓我成為無根的浮萍,談慶明老師建議鄭先生讓我跟一個研究圓柱尾流混沌現象的實驗團隊一起做研究.這個項目的立項,白老師和曹師兄起了很大作用.團隊由林貞彬和鄂學全兩位老師領銜,成員中有位動手能力極強的年輕人,李東輝.團隊研究的目的是通過測量和分析圓柱尾流固定點速度分量的時間序列來厘清圓柱尾流的初始轉捩區域是否存在混沌現象,若是,再繼續研究混沌在湍流發展過程中的作用.

自1963 年Lorenz[8]發表Lorenz 混沌方程至20 世紀八九十年代,混沌現象的基本特性及由規則運動向混沌運動轉化的機制都已清楚.前者包括分數維特征和系統對初值的敏感性,后者包括周期倍化至混沌[9]、準周期至混沌[10]及間歇性至混沌[11]等路徑.這里只粗略介紹一下最基本的概念和現象,更詳細的請見文獻[1,12].系統對初值的敏感性指的是小擾動隨時間呈指數增長.令d(0)為0 時兩條任意軌跡之間的小距離,d(t)為時間t時它們之間的距離,對于真正的低維確定性混沌,我們有

其中,λ1被稱作最大的Lyapunov 指數.混沌吸引子的軌跡在相空間中有界.相鄰軌跡的指數發散導致的不斷拉伸,加上吸引子的有界性引起的不時折疊,使混沌吸引子經常成為分形,其特征是

其中,N(ε) 表示在相空間中完全覆蓋吸引子所需的線性長度不大于 ε 的(最小)盒子數.D稱為吸引子的計盒維數.通常它是一個非整數.對于混沌 Lorenz吸引子來說,D為2.06.值得指出,方程(2)也適用于度量無人經營過的山路或海岸線的長度.這種情況,D通常大于1 但小于2.容易看出,當 ε 越來越小時,山路或海岸線的長度越來越長.特別地,當 ε →0 時,山路或海岸線的長度也變成無限長.

當時,混沌時間序列分析的基本框架也已完善.給定一個時間序列x(1),x(2),···,x(n),首先需構造一個相空間,其中的向量為

其中,m稱作嵌入維,L是延遲時間.然后,可以用Wolf 等[13]的算法計算最大Lyapunov 指數,并用Procaccia 和Grassberger[14]的算法計算關聯維(計盒維數D的一個很緊湊的下限)和關聯熵(λ1的一個很好估計).一般假設,當關聯維為非整數、Lyapunov 指數和關聯熵為正時,所研究的時間序列是混沌的.當我們研究圓柱尾流有無混沌現象時,也用到了這個假設[15].但是這樣做是不對的.有諸多原因.其一是,用Wolf 等[13]的算法計算最大Lyapunov指數時,小擾動的指數增長特性,只是被假設,但沒被證明.事實上可以證明,從任意一組隨機的數據,都可以計算出一個正的Lyapunov 指數.一方面,這與混沌來自決定性系統這個原理完全不一致.另一方面,這并不表示隨機系統里小擾動有指數發散的特征.事實上,證明一個復雜時間序列的小擾動有指數發散的特征,是一個極富挑戰性的問題.第二個原因是,一類很常見的隨機過程,1/fα過程,也擁有非整數的關聯維和正的關聯熵,因而會被誤判成混沌[16-17].第三個重大原因是,實驗數據總有噪聲.當時的分析方法沒法幫助研究者撇開噪聲,把數據里最重要的特性刻畫出來.第四個原因是技術上而非原理性的:實驗數據長度有限,因此式(3) 里m和L的選擇就變得很重要.這個問題被稱作相空間的最優重構.1993 至1994 年間,我和鄭先生發表了三篇文章,把這些問題都解決了.下面解釋一下基本原理.我們從相空間的最優重構講起.

從原理上來說,式(3) 里的m只要夠大就行,L可以任意.m夠大指的是不小于2 倍的計盒維數D[18].問題在于D需要估計而非已知.所以相空間的最優重構這個問題,本質上是在數據有限這個約束條件下,估計最小可用的嵌入維m,然后選擇L使得動力學的演化能被最清楚地觀察到.m不能太小,可以通過考慮將南北極投影到赤道來理解-在2 維的赤道上,沒法區分南北極,所以重構的相空間必須夠大,使得地球看起來確實像個地球.L的選擇與相空間里運動的速率有關-若速率變化忽大忽小非常極端,則動力學的演化很難被觀察清楚.當時已有一個很好的方法,假的最近鄰(false nearest neighbor)方法,來最優地重構相空間[19-20].這是一個幾何的方法,其基本想法是,當m增大至最小可用的嵌入維時,假的最近鄰的個數會急劇減少.嵌入維選定后,延遲時間可以通過讓假的最近鄰個數最小而得到.

我和鄭先生構造的方法是動態的[21-23].基本想法是假定最近鄰通過動態迭代后的發散遠快于真的最近鄰的發散.用方程表述是

其中,Vi和Vj等是重構的向量,k是演化時間,尖括號表示滿足下述條件的所有可能的 (Vi,Vj) 對的系綜平均

其中,εk和 Δ εk是任意選擇的小距離.這里的每一個不等式對應一個球殼,引入球殼的根本作用是為了刻畫噪聲的影響.這一點在后面討論混沌的直接動力學判據時,會更清楚.為最優地決定m和L,可以固定k至一個較小的值.以Rossler 吸引子的重構為例,見圖1,我們發現,m=3,L=8 是一個優化解.更多的細節請見文獻[21].

圖1 Rossler 吸引子的重構Fig.1 Reconstruction of Rossler's attractor

值得指出,這個動態的方法和假的最近鄰法,當時是,現在仍是解決相空間重構的兩個最系統的方法.我們的動態方法提供更多的信息.以Lorenz 吸引子Lyapunov 指數的估計為例,見圖2.當向量對(Vi,Vj)無所限制,即下面不等式(6)里的w為1 時,Λ(k)隨時間變化的斜率比Lyapunov 指數小.原因是,對應Lyapunov 指數為0 的軌道上的切向運動也被計算在里邊了.讓w大一些,如選取54,則 Λ (k) 隨時間的變化的斜率與最大Lyapunov 指數完全一致

圖2 Lorenz 吸引子Lyapunov 指數的估計Fig.2 Estimation of the Lyapunov exponent for the Lorenz attractor

更重要的是,這個動態的方法給出了一個混沌的非?？煽康闹苯觿恿W判據.如圖3(a)所示,對真正的混沌系統來說,來自不同球殼的隨時間變化的指數曲線會形成一個包絡線,其斜率給出最大Lyapunov 指數的一個非常好的估計.對隨機數來說,如圖3(b)所示,來自不同球殼的隨時間變化的指數曲線,雖然都能給出一個正的Lyapunov 指數(相當于L(k)/k),但因為是分散的,用不同的球殼去估計最大Lyapunov 指數時,得到的結果都不一樣.這是隨機系統的本質特征.特別地,當數據集的大小趨向于無限,而球殼的大小趨向于0 時,最大Lyapunov 指數將變得無限大.這也是隨機系統的一個根本特征——熵無限大的一個很好的體現.

圖3 隨時間變化的指數曲線作為混沌的直接動力學判據:(a) Lorenz 系統,(b) 隨機系統.數字從小到大對應球殼從大到小[21]Fig.3 Time-dependent exponent curves for the chaotic:(a) Lorenz data and (b) IID random variables,where the curves,from bottom up,correspond to shells from larger to smaller[21]

有了這個混沌的直接動力學判據,我們可以搞清楚圓柱近尾流低雷諾數的流動是否有混沌運動[24].圖4(a)顯示的是Re=136 時的相圖.類似的相圖在文獻里經常出現,如有節律的化學反應和原子力顯微鏡[25],且被當作混沌運動存在的證據.這里,小擾動以冪律而非指數方式增長,l nεt～t1/2,見圖4(b),所以這是極限環上的布朗運動而不是混沌運動.

圖4 圓柱近尾流低雷諾數(Re=136)流動的混沌行為:(a) 相圖,(b)lnt=ln||Vi+k-Vj+k||隨時間的變化[24]Fig.4 Test of chaotic behavior in the near wake cylinder of low Reynolds number (Re=136):(a) phase diagram and (b) temporal evolution of l nt=ln||Vi+k-Vj+k|| [24]

懂了混沌的直接動力學判據后,就能很方便地理解噪聲對動力系統的影響.對混沌系統來說,噪聲的作用是破壞圖3(a)所示的包絡線,隨著噪聲越來越大,包絡線從上到下,對應球殼從小到大,將逐漸被破壞,直至包絡線完全消失,系統小擾動的指數發散特征不再能被分辨出來.這個現象可以用來估計噪聲的大小[26].總的來說,這套方法能非常好地刻畫噪聲對動力系統,包括半導體激光系統的動力學的影響[27],也能非常好地研究一個相反的問題——噪聲引發的混沌[28-30].細節就不再贅述了.這里只強調一下,這也是我在上一節提到的,能在沒有任何基礎的情況下,簡單融入劉教授的團隊,研究用半導體激光做混沌安全通訊的原因.

3 基于依賴指數的李雅普諾夫指數(SDLE)的多尺度分析

混沌理論對20 世紀科學思想的最大沖擊是,隨機性可以不由無限維的隨機系統產生,而由有限維的確定性系統產生.這個將隨機系統與確定性系統對立的表述,產生的一個很嚴重的后果是,早期的一些研究者,包括我自己,在選擇了混沌研究后,會忽略隨機系統的研究.1997 年在UCLA 的電子工程系研究互聯網上的交通流時,我意識到這個偏見很糟糕.任一時間、任一地點的互聯網牽涉到的用戶數,都遠多于地面交通牽涉到的車輛數.所以正常情況下互聯網上的交通流必然是隨機的,這種情況是相對于網絡受攻擊時的狀態而言,當網絡受攻擊時,其上的交通流會有一個成分,由程序(即規則)產生,因此會不再是完全隨機的.互聯網上交通流的一個基本特征是長程相關性[31].雖然離散的混沌映射能模擬這個特征[32],但兩個著名的隨機分形模型,1/fα過程和乘法級聯多分形模型,能更全面地模擬互聯網上交通流的特征[33-36].值得注意,乘法級聯多分形模型是為理解湍流而被發展起來的[37-41],且至今仍是湍流最好的模型之一.早期復雜性科學研究的一大目標是用低維混沌理論研究湍流.所以這里一個深刻的問題是如何將低維混沌與隨機分形理論有機地結合起來.SDLE 是為了實現這個目標的一個很有成效的嘗試.為了充分理解SDLE 的涵義,下面先介紹一下 1/fα過程.隨機系統的其他主要模型,包括ARIMA 模型,Levy 過程,和乘法級聯多分形模型,就不在此贅述了(可見文獻[1]).

在表征復雜系統的活動類型中,一類非常普遍且令人費解的行為是 1/fα過程,其功率譜密度隨頻率(時域)或波長(空間域)以冪律方式衰減[1],因此其維數不能通過主成分分析等方法來有效降低[42].這種過程的一個子類,通常表示為 1/f2H+1過程,具有長程相關性,其中,H稱作赫斯特指數(Hurst Parameter).取決于 0＜ H ＜1/2,H=1/2,或 1/2 ＜ H ＜1,它們分別被稱作具有反持久相關性、無記憶或僅短程相關性,及持久長程相關性(長記憶)[43].此類過程的突出例子包括視覺[44]、金融[45-46]、DNA 序列[47-51]、人類認知[52]、全球恐怖主義[53]、協調[54]、姿勢[55-56]、心臟動力學[57-61]以及素數的分布[62]等.

具有長記憶過程的精確定義如下:均值為 μ、方差為 σ2,且協方差平穩的隨機過程X={X(t),t=0,1,2,···},當其自相關函數r(w),w≥0 以冪律方式衰減時[63]

則有長程相關性,其中 0 ＜H＜ 1.當 1/2 ＜H＜1,這是為什么這種過程被稱作具有長記憶.過程X的功率譜密度 (PSD) 為它的積分被稱為隨機游走過程,其PSD 為作為1/f過程,它們不能被馬爾可夫過程或ARIMA 模型恰當地建模[64],因為這些過程的 PSD 與1/f明顯不同.為了充分模擬1/f過程,必須使用分數階過程.最流行的模型是分數布朗運動模型[43].

值得強調的是,充分發展的湍流的Kolmogorov 能譜對應H=1/3,其效應在用雷達觀察到的海平面的海雜波上也有反映[65-66].

在上述基礎,我們可以討論SDLE.

假設一個適當的相空間已經重構好了,我們考慮一個軌跡的集合.將兩個相鄰軌跡之間的初始分離表示為 ε0,它們在時刻t和t+Δt的平均距離為εt和 εt+Δt.相空間中軌跡分離的一個示意如圖5 所示.我們考察 εt和 εt+Δt的關系.當 Δt→0,我們有

圖5 相空間中軌跡分離的示意圖Fig.5 Schematic of evolution of separation of trajectories in phase space

其中,λ (εt) 就是 SDLE,由下式給出

也可以把它表示成一個微分方程

有趣的是,λ (εt) 的計算,完全可以借用上一節描述的我和鄭先生發展的方法.具體是,還是基于不等式(5)描述的球殼,我們直接得到

其中,t和 Δt是采樣時間的整數倍,尖括號表示球殼內所有下標為i,j的平均值.注意,SDLE 蘊含Lyapunov 指數的概念,但比后者的概念寬廣,因SDLE 是一個函數,而Lyapunov 指數只是一個數.

粗略地說,當時我和鄭先生發展的方法可以稱作是一個積分的形式,SDLE 的表述的是一個微分的形式.由積分形式轉至微分形式帶來的好處是,在SDLE 的框架下,對所有已知的時間序列模型,我們都可以證明它們有其獨特的標度律.下面列舉一些,更多的細節請見文獻[1,61-67].

性質1:對確定性混沌,可以證明,在小尺度上

性質2:對于有噪聲的混沌和噪聲引起的混沌,在小尺度上,有

其中 γ ＞0 決定信息衰減的速率;

性質3:對1/f2H+1過程,可以證明

性質4:對 α-stable Levy 過程,可以證明

性質5:對隨機振動來說,依賴于相空間的重構,λ(ε)≈-γlnε 和 λ (ε)≈-Hε-1/H都可被觀察到;

性質6:對于具有多尺度行為的復雜運動,上述不同的標度律有可能在不同的 ε 范圍內都被觀察到.

上述6 個性質,涵蓋了幾乎所有已知的時間序列模型.因而我們可以期望,SDLE 能統一所有已知的復雜性的度量.事實上也確實如此.比如在刻畫腦電動力學時,我們發現,常用的復雜性的度量,包括Lyapunov 指數,關聯維,關聯熵,赫斯特指數,等都可以被對應到SDLE 在某個尺度的值.細節請見文獻[1](15.7.1 節)和文獻[68].因這個原因,SDLE 在臨床腦電的分析中有極大的應用前景[69].

因性質6 有重大意義,我們用一個例子使其具體化.考察一個動力系統

其中,[xn] 表示xn的整數部分,ηn是在區間 [-1,1] 上均勻分布的隨機數,σ 刻畫噪聲的強度,F(y) 是個由下述方程給出的函數

函數F(y) 描述的是一個混沌動力學,其Lyapunov指數為 l n(2+Δ).所以,由方程(16)描述的動力系統是一個整數格點上的隨機運動加上整數格點內的混沌運動.因此,當沒有噪聲時,我們應該在小尺度上觀察到性質1,而在較大尺度上觀察到性質3.確實如此,如圖6 所示.當有噪聲時,性質1 被性質2 替代.這里,圖6(a)用了5000 點,圖6(b)用了500 點.后一種情形,因點太少,整數格點上的布朗運動沒能被刻畫好(事實上,真正牽涉整數格點的個數只有幾十).總的來說,用這么少的點刻畫這么豐富的動力學,很令人驚訝.

圖6 小尺度上的混沌運動與大尺度上的布朗運動Fig.6 Chaotic motion on small scales and Brownian motion on large scales

性質6 的一個直接后果是,我們可以把一個混沌或分形時間序列拆成很多段,中間插入周期性運動,則原來的混沌或分形特征還是能被SDLE 很好地刻畫出來.事實上,前者對應間歇性混沌.一個例子由圖7 給出.分形+振蕩是長時間監測的心跳動力學的一個基本特征.這里,振蕩由呼吸引起.這樣的分形很難被SDLE 以外的任何方法刻畫.細節就不再此贅述了,請見文獻[61].

圖7 (a) 由Logistic 映射 xn+1=axn (1-xn) 生成的間歇性混沌時間序列,a=3.8284.(b) 由10000個點的時間序列計算出的 SDLE 曲線,m=4,L=1,最大的球殼尺寸為 (2 -13.5,2 -13).由箭頭 A 指示的混沌標度區清晰可見.箭頭 B 和 C 表示的區域對應周期性運動和混沌運動之間的過渡,箭頭 D 表示大尺度的振蕩運動Fig.7 (a) Intermittent chaos generated by Logistic map xn+1=axn (1-xn),a=3.8284.(b) SDLE curve calculated from a time series of 10000 points,m=4,L=1,and the largest shell is (2 -13.5,2 -13).The chaotic scale region indicated by arrow A is clearly visible.The regions represented by arrows B and C correspond to the transition between periodic and chaotic motions,while arrow D represents large-scale oscillatory motion

思考一下,我們就能夠明白,SDLE 的提出,拓展了傳統混沌動力學的研究范疇,因后者一般只研究由不變測度刻畫的含有混沌吸引子的動力系統.在大數據時代,數據一般都是被連續觀察的.真正讓人激動的數據是那些含有狀態突變的數據.這種情形,雖不易由大多數已知的方法來處理,但對SDLE 并沒造成任何挑戰.這里,值得提一下河流徑流動力學.許多研究者猜測河流動力學是混沌的,但分析結果不能讓人信服.原因是,河流動力學是強間隙性的—旱季的枯水期和雨季的豐水期,其動力學是完全不一樣的.因此,河流的流動形態中若存在混沌,該混沌必須是間隙性的.所幸,這里的混沌很容易通過SDLE 刻畫.一個例子如圖8 所示.

圖8 美國Colorado 河的間隙性混沌.這里,紅色的曲線代表由下節的非線性自適應濾波器處理過的數據Fig.8 Intermittent chaos in the Colorado River,USA.Here,the red curves represent the data processed by the nonlinear adaptive filter to be explained in the next section

4 自適應去趨勢、去噪、多尺度分解和分形分析

2006 年,美國因在伊拉克和阿富汗缺兵力,由陸軍研究部負責,請一家公司設計了一款能用于野外很方便采集腦電的頭盔,然后請杜克大學醫學院的一個團隊采集了一組腦電,最后找到我們,希望我們能設計出一套通過腦電分析,評估受炮彈沖擊波沖擊而腦震蕩的士兵,是否已恢復可以再上前線了.因這個經費來之不易,我們做的非常認真.做了幾個月,發現實驗得到的腦電信號,其實沒有任何真正腦電的信息.但我們不能就這樣告訴陸軍研究部.所以我們設計了一個新的分析方法,證明給我們的腦電,只有頭部搖晃的信號,加上電源60 Hz 的信號(中國是50 Hz),及一些噪音.陸軍研究部雖然對我們的分析很滿意,但畢竟很失望,因為設計的頭盔完全無效.我自己得到的一個深刻體會是,實驗數據采集的團隊,必須包含能做分析的人員,不然采集的數據可能都是垃圾.

我們的方法是這樣的[70-75].先把一個時間序列分成長度為奇數w=2n+1 的小段,相鄰兩段重疊n+1個點.這樣做引入了一個時間尺度=(n+1)τ,其中 τ 為采用時間間隔.對于每一段,我們擬合一個M階的最佳多項式.請注意,M=0 和 1 分別對應于分段常數和分段線性擬合.分別用yi(l1) 和yi+1(l2)表示擬合的第i段和第(i+1)段的多項式,其中l1,l2=1,2,···,2n+1.注意最后一段的長度可能小于2n+1.我們將重疊區域的擬合定義為

其中,w1=1-(l-1)/n,w2=(l-1)/n可以重寫為1-dj/n,j=1,2,其中dj表示那一點與y(i) 和y(i+1)的中心的距離.這意味著權重隨點和線段中心之間的距離線性減小.這種加權確保了對稱性并有效地消除了相鄰段邊界周圍的任何跳躍或不連續性.實際上,該方案確保擬合處處連續,在非邊界點處光滑,在邊界處具有右導數和左導數.該方法能有效找出任何一種信號的趨勢,假如存在的話.

回到杜克大學的腦電數據.兩個例子如圖9 所示.紅色的曲線由我們的自適應濾波器得到.它們對應頭部的晃動.較細的黑色趨勢線由另一個很有名的濾波器Loess[76]得到.Loess 等價于Savitzky-Golay 濾波器[77].可以看出,紅色的趨勢線比黑色的更精確.

圖9 腦電因頭部晃動引起的趨勢Fig.9 Trend of EEG caused by head shaking

值得比較這個濾波器與前面提到EMD[2].圖10(a)顯示的是全球年海表溫度 (SST)及線性、拋物線、非線性擬合,右邊顯示的是相應的殘差序列.圖10(b)與文獻[2]的結果完全一致.值得注意,EMD 是個二分法,時間尺度每次減半.這里的自適應算法牽涉到的時間尺度是連續的,因此,在許多應用場合,自適應算法應該比EMD 更準確.

圖10 (a)全球年海表溫度(SST)及線性、拋物線、非線性擬合和(b)相應的殘差序列Fig.10 (a) Global annual sea surface temperature (SST) and linear,parabolic and nonlinear fitting and (b) corresponding residual series

自適應濾波器有非常好的去噪功能.一個例子如圖11 所示.可以看出,自適應濾波器比基于混沌的方法和小波變換都好.細節請見文獻[71].

圖11 混沌洛倫茲信號的去噪Fig.11 Denoising of a chaotic Lorenz signal

自適應濾波器也有極好的信號提取功能.一個例子如圖12 所示.其中,圖12(a)和12(b) 是以采樣頻率為250 Hz 采集的 EEG 和 ECG 信號;圖12(c)和12(d)是通過使用自適應和小波變換從 EEG 中提取的ECG 成分.圖12(a)中的箭頭突出顯示了 EEG信號中包含的 ECG 分量.在呼吸綜合征的監測和治療中,腦電EEG 的分析起到很重要的作用.但混雜了心電ECG 成分的腦電,使得通過分析判斷是否有呼吸綜合征變得很困難.因此,一個重要步驟是先去掉心電ECG 的成分.如圖所示,自適應方法比小波變換更優.

圖12 (a,b)以采樣頻率為250 Hz 采集的 EEG 和 ECG 信號,單位為mV;(c,d)使用自適應和小波變換從 EEG 中提取的ECG 成分.(a)中的箭頭突出顯示 EEG 信號中包含的 ECG 分量Fig.12 (a,b) EEG and ECG signals acquired at a sampling frequency of 250 Hz in mV;(c,d) ECG components extracted from EEG using adaptive and wavelet transforms.Arrows in (a) highlight ECG components included in the EEG signal

理解了自適應濾波器上述功能后,就很容易理解它的多尺度分解的功能.這相當于用一系列窗口w1,w2,...,wk,得到一連串趨勢線,任意兩個趨勢線的差,都代表原來時間序列在對應那兩個窗口的頻段的成分.對具體例子感興趣的,請見文獻[73].

最后,我們解釋如何用自適應濾波器作分形和多分形分析.事實上,當數據有趨勢時,它是所有估計赫斯特指數及廣義的赫斯特指數譜的方法中最好的.這里只給出一些最基本的結果.細節請見文獻[73].

然后,我們用窗口大小為w對u(n) 擬合一個全局趨勢v(n).殘差u(i)-v(i) 表征圍繞全局趨勢的波動,其方差是赫斯特指數H的一個非常好的估計

將上式的平方與開根號稍作變化,我們就得到估計廣義的赫斯特指數譜的多分形分析

將式(20)應用于兩個隨機游走過程u(n) 和z(n),我們可以計算它們間互相關的長程相關指數Huz

圖13(a)顯示的是一個典型的腦電信號的分形分析結果.從中我們可以定義三個參數,短時間和長時間的赫斯特指數(記為HS和HL),及飽和度(saturation level),就可以把三類腦電數據,正常(Group H),癲癇發作(Group S),癲癇病人未發作(Group E),幾乎完全區分開.近期,我們發現這個方法有很大的臨床應用前景[78].

圖13 (a)一個典型腦電信號的分形分析結果;(b) 三類數據,正常(Group H)、癲癇發作(Group S)、癲癇病人未發作(Group E)的分類,精度幾乎100%Fig.13 (a) Typical adaptive fractal analysis of a an EEG signal;(b)Classification of three types of EEG data,normal (Group H),seizures(Group S),and seizure patients not having seizure (Group E),with an accuracy of almost 100%

5 結語

5 G 時代的到來及 6 G 時代的快速臨近,科學、工程和社會中積累的大數據將很快變得更加巨大.任何人都不得不抓住這個前所未有的機會.雖然計算機科學家們正在努力開發更強大的數據庫管理和機器學習方法,但現在是進入下一階段的時候了.這個階段將以捕獲大數據背后的動力學及深刻理解人腦的工作機制為基本目的.就數據分析而言,我們可以洞見,機器學習和復雜性科學的主流方法在未來會相互補充并互相促進.為了幫助加速這個結合,我們提倡協同使用基于主流機器學習的方法和復雜性科學中的多尺度方法.

最后我們強調,在科學研究中,最重要的是一直保有一顆好奇、不斷探索的心.只有這樣,才能發現真正重大的科學問題并切實解決實際問題.這種普世的精神,鄭先生等老一輩科學家們一直在言傳身教.望讀者們能和我們共勉！