簡(jiǎn)介這份MATLAB源碼包面向信號(hào)處理、非線性動(dòng)力學(xué)與時(shí)間序列分析研究者圍繞三維相空間重構(gòu)PSR提供一套完整可運(yùn)行的算法實(shí)現(xiàn)適用于科研復(fù)現(xiàn)、教學(xué)演示與工程項(xiàng)目開(kāi)發(fā)。壓縮包共11個(gè)文件以4個(gè)M腳本為核心涵蓋互信息法延遲估計(jì)、FNN嵌入維計(jì)算、Lorenz混沌系統(tǒng)時(shí)間序列生成與三維可視化另含3張結(jié)果圖片、1份txt數(shù)據(jù)、1個(gè)Markdown說(shuō)明文檔與License授權(quán)文件包體僅205KB。目前已有183人學(xué)習(xí)下載。通過(guò)這套代碼可直觀理解Takens延時(shí)嵌入定理的實(shí)操流程從原始一維時(shí)間序列到重構(gòu)相空間再到吸引子形態(tài)繪制完整呈現(xiàn)相空間重構(gòu)的關(guān)鍵步驟代碼結(jié)構(gòu)清晰參數(shù)選擇部分便于替換到真實(shí)信號(hào)中進(jìn)一步探索是學(xué)習(xí)混沌分析的良好參考實(shí)現(xiàn)。1. 一份 PSR 三維重構(gòu)源碼能幫你看到什么“所有代碼_psr_三維重構(gòu)_相空間_相空間重構(gòu)_straightxx8_源碼”是典型的資源站下載包一堆關(guān)鍵詞拼起來(lái)的壓縮包沒(méi)有說(shuō)明書也沒(méi)有版本號(hào)。它的技術(shù)主線很聚焦——用相空間重構(gòu)Phase Space ReconstructionPSR把一維時(shí)間序列映射到三維空間里把混沌信號(hào)藏在時(shí)間軸里的吸引子結(jié)構(gòu)“畫”出來(lái)。這類源碼在振動(dòng)故障診斷、生理信號(hào)分析、非線性時(shí)間序列預(yù)測(cè)里出現(xiàn)頻率很高適合手里有一段實(shí)測(cè)數(shù)據(jù)、想判斷它到底是隨機(jī)噪聲還是確定性結(jié)構(gòu)或者想給分類模型造一個(gè)更好分特征的人。它解決的核心問(wèn)題可以用一句話概括同一段波形時(shí)域里看不出規(guī)律升到三維空間后規(guī)則結(jié)構(gòu)立刻現(xiàn)形。下面我按復(fù)現(xiàn)這類源碼包的順序把原理、算法、實(shí)現(xiàn)、踩坑和定量分析一次講清。2. 相空間重構(gòu)原理與參數(shù)選型τ 和 m 為什么決定三維圖長(zhǎng)什么樣相空間重構(gòu)在大部分人聽(tīng)來(lái)像玄學(xué)核心其實(shí)是一句話一段標(biāo)量時(shí)間序列里藏著系統(tǒng)全部狀態(tài)的演化軌跡。決定三維重構(gòu)圖好不好看的只有兩個(gè)參數(shù)——延遲時(shí)間 τ 和嵌入維數(shù) m。源碼包里幾乎所有子程序都在圍著這兩個(gè)參數(shù)轉(zhuǎn)讀懂了它倆任何 PSR 源碼都不會(huì)再看暈。2.1 Takens 嵌入定理從一維序列恢復(fù)吸引子拓?fù)銽akens 在 1981 年證明的嵌入定理是這套方法的根基。假設(shè)原始動(dòng)力系統(tǒng)是 d 維的我們能觀測(cè)到的只是其中一個(gè)坐標(biāo)的采樣序列 x(t)。構(gòu)造延遲向量X(t) [x(t), x(tτ), ..., x(t(m-1)τ)]當(dāng)嵌入維數(shù) m ≥ 2d1 時(shí)重構(gòu)后的軌跡與原始吸引子是拓?fù)涞葍r(jià)的。翻譯成人話雖然每個(gè)時(shí)刻只觀測(cè)到一個(gè)數(shù)值但把“現(xiàn)在”和“未來(lái)幾個(gè)時(shí)刻”拼成一個(gè)向量足夠還原系統(tǒng)內(nèi)部狀態(tài)的演化關(guān)系。所謂三維重構(gòu)就是取 m3 的特例三個(gè)坐標(biāo)軸分別是 x(t)、x(tτ)、x(t2τ)。這里有個(gè)特別容易被誤解的點(diǎn)重構(gòu)坐標(biāo)沒(méi)有物理單位也不代表原系統(tǒng)里的某個(gè)物理量它只是延遲副本構(gòu)成的抽象空間。所以別給坐標(biāo)軸硬標(biāo)“電壓”“位移”之類的量。拓?fù)涞葍r(jià)的意義在于幾何不變量可以保留——關(guān)聯(lián)維數(shù)、Lyapunov 指數(shù)這些反映系統(tǒng)本質(zhì)的量在重構(gòu)空間里算和在原系統(tǒng)里算結(jié)果一致這是后面做定量分析的前提。工程上d 一般未知所以 m 通常從 2 試到 10 左右看結(jié)構(gòu)和指標(biāo)是否穩(wěn)定。如果只是想“畫個(gè)三維吸引子看看”m 固定為 3 就夠了。源碼包里大量出現(xiàn) m3 不是偷懶是可視化場(chǎng)景下的合理選擇。2.2 延遲時(shí)間 τ 的兩種算法自相關(guān)法與互信息法τ 選小三個(gè)坐標(biāo)高度相關(guān)軌跡擠成一條線τ 選大三個(gè)坐標(biāo)近似獨(dú)立軌跡變成隨機(jī)點(diǎn)云。自相關(guān)法和互信息法是最常用到的兩種選法。自相關(guān)法算的是 x(t) 和 x(tτ) 之間的線性相關(guān)系數(shù)隨 τ 的衰減常見(jiàn)準(zhǔn)則取第一次降到 1/e 的位置作為 τ。優(yōu)點(diǎn)是快缺點(diǎn)是只捕捉線性依賴對(duì)非線性結(jié)構(gòu)不敏感。import numpy as np def autocorr_tau(signal, stop1.0 / np.e): x signal - signal.mean() n len(x) # 補(bǔ)零到 2n用 FFT 算線性自相關(guān)避免逐點(diǎn)循環(huán) fft_x np.fft.fft(x, n2 * n) acov np.fft.ifft(fft_x * np.conj(fft_x)).real[:n] / n acov acov / acov[0] # 歸一化到 lag0 時(shí)相關(guān)系數(shù)為 1 for tau in range(1, n): if acov[tau] stop: return tau return n - 1邏輯說(shuō)明先減均值消除直流分量再補(bǔ)零做 FFT 計(jì)算自相關(guān)比逐點(diǎn)雙重循環(huán)快幾個(gè)數(shù)量級(jí)。除以 acov[0] 完成歸一化閾值就直接用 1/e。對(duì) Lorenz 這類信號(hào)在 dt0.02 時(shí)算出的 τ 通常在 5 到 15 之間和文獻(xiàn)里“延遲時(shí)間取自相關(guān)第一次過(guò)零點(diǎn)附近偏小一點(diǎn)”的經(jīng)驗(yàn)吻合。保守的寫法是取第一次過(guò)零點(diǎn)但那給出來(lái)的 τ 往往偏大軌跡會(huì)明顯變稀疏。自相關(guān)法的局限在于它只度量線性相關(guān)性?；バ畔⒎▌t能捕捉非線性依賴它把信號(hào)值域分成若干格子統(tǒng)計(jì)滯后 τ 的兩個(gè)變量共享多少信息量取第一極小點(diǎn)作為 τ。def mutual_information(signal, tau, bins32): x signal[:-tau] y signal[tau:] lo, hi np.min(signal), np.max(signal) # 聯(lián)合直方圖固定使用全序列值域保證不同 tau 之間可比 cxy, _, _ np.histogram2d(x, y, binsbins, range[[lo, hi], [lo, hi]]) n cxy.sum() pxy cxy / n px pxy.sum(axis1) py pxy.sum(axis0) mi_val 0.0 for i in range(bins): for j in range(bins): if pxy[i, j] 0: mi_val pxy[i, j] * np.log(pxy[i, j] / (px[i] * py[j])) return mi_val def mi_first_min(signal, tau_max80, bins32): vals [mutual_information(signal, t, binsbins) for t in range(1, tau_max 1)] for i in range(1, len(vals) - 1): if vals[i] vals[i - 1] and vals[i] vals[i 1]: return i 1 # 索引 i 對(duì)應(yīng) tau i1 return tau_max參數(shù)說(shuō)明bins 取 32 是常見(jiàn)折中數(shù)據(jù)總量少于幾千點(diǎn)時(shí)降到 16否則聯(lián)合直方圖大量格子為零互信息抖動(dòng)很厲害。tau_max 要覆蓋信號(hào)的一個(gè)主周期dt0.02 的 Lorenz 軌道時(shí)間常數(shù)在 1 秒量級(jí)tau_max 取 80 足夠。代碼里返回的是第一個(gè)局部極小點(diǎn)不是全局最小點(diǎn)這是 Fraser-Swinney 方法的經(jīng)典約定。兩套算法結(jié)果不一致時(shí)怎么辦比如自相關(guān)給 8、互信息給 15先畫互信息曲線看第一極小是否明顯再在兩值之間取偏大者做可視化。稍大的 τ 能把軌跡拉開(kāi)、看到更多折疊結(jié)構(gòu)如果差異超過(guò) 3 倍多半是信號(hào)有趨勢(shì)或周期性太強(qiáng)先去趨勢(shì)再說(shuō)。2.3 嵌入維數(shù) m 的確定從偽近鄰到“夠用就好”如果只是三維可視化這部分可以跳過(guò)。但源碼包通常還帶 G-P 算法或偽近鄰法說(shuō)明作者意圖不止畫圖。偽近鄰的思路在 m 維空間里一個(gè)點(diǎn)的大部分近鄰應(yīng)該是“真鄰居”如果升到 m1 維后原本的近鄰跑遠(yuǎn)了說(shuō)明那些是低維投影造成的假鄰居。m 從 1 遞增偽近鄰比例降到接近 0 時(shí)的 m 就是合適嵌入維。G-P 算法從另一個(gè)方向逼近在重構(gòu)空間里統(tǒng)計(jì)距離小于 r 的點(diǎn)對(duì)比例得到關(guān)聯(lián)積分 C(r)log-log 坐標(biāo)下無(wú)標(biāo)度區(qū)的斜率就是關(guān)聯(lián)維數(shù) D2。隨著 m 增大確定性混沌系統(tǒng)的 D2 會(huì)飽和在某個(gè)值附近如果 D2 一直漲信號(hào)大概率是隨機(jī)噪聲。僅這一條就常被用來(lái)區(qū)分“混沌”和“純隨機(jī)”。工程選型建議畫三維圖用 m3估算關(guān)聯(lián)維數(shù)或 Lyapunov 指數(shù)用 m5 到 7 起步高維系統(tǒng)通常要 m≥8。m 不是越大越好——樣本量固定時(shí)空間維數(shù)越高數(shù)據(jù)越稀薄距離估計(jì)全部失真。經(jīng)驗(yàn)上要求重構(gòu)后的點(diǎn)數(shù) n-(m-1)τ 至少是 m 的 10 倍否則后面算關(guān)聯(lián)維數(shù)基本是噪聲。3. 從源碼包到最小復(fù)現(xiàn)Lorenz 信號(hào)的三維相空間重構(gòu)拿到這類源碼包最常見(jiàn)的做法是先把環(huán)境配干凈用一段已知答案的混沌信號(hào)把流程跑通再換自己的數(shù)據(jù)。不要一上來(lái)就上真實(shí)信號(hào)因?yàn)檎鎸?shí)信號(hào)里的噪聲和趨勢(shì)會(huì)讓“圖不對(duì)”時(shí)無(wú)法判斷是自己錯(cuò)了還是數(shù)據(jù)本身有問(wèn)題。3.1 造一段已知答案的測(cè)試信號(hào)Lorenz 系統(tǒng)from scipy.integrate import solve_ivp def lorenz(t, state, sigma10.0, rho28.0, beta8.0 / 3.0): x, y, z state return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z] dt 0.02 t_end 120 t_eval np.arange(0, t_end, dt) sol solve_ivp(lorenz, [0, t_end], [1.0, 1.0, 1.0], t_evalt_eval, methodRK45, rtol1e-8) x sol.y[0] x x[2000:] # 丟掉前 40 秒瞬態(tài) print(f剩余點(diǎn)數(shù): {len(x)})邏輯說(shuō)明Lorenz 方程在 sigma10、rho28、beta8/3 的經(jīng)典參數(shù)下處于蝴蝶混沌區(qū)初值隨便給只要不落在平衡點(diǎn)附近就行。積分完成后把前 2000 點(diǎn)丟棄因?yàn)閺某踔碉w到吸引子上的過(guò)渡段會(huì)在重構(gòu)圖里多出一條“飛線”。兩個(gè)參數(shù)要記住dt 是采樣間隔直接決定 τ 的物理含義rtol1e-8 防止數(shù)值誤差讓軌跡跳到另一個(gè)分支。真實(shí)數(shù)據(jù)沒(méi)有積分這一步但一定有采樣率建議一開(kāi)始就把 τ 的離散值換算成物理時(shí)間。提示真實(shí)信號(hào)做相空間重構(gòu)前先確認(rèn)采樣率和主頻帶。采樣率過(guò)高時(shí)先降采樣否則重構(gòu)點(diǎn)數(shù)暴漲圖也卡τ 的物理意義也容易算錯(cuò)。3.2 相空間重構(gòu)核心實(shí)現(xiàn)與三維可視化def psr_reconstruct(signal, tau, m3): n len(signal) rows n - (m - 1) * tau if rows 0: raise ValueError(n-(m-1)*tau 為負(fù)數(shù)據(jù)太短或參數(shù)太大) mat np.empty((rows, m)) for i in range(m): mat[:, i] signal[i * tau : i * tau rows] return mat tau 12 mat psr_reconstruct(x, tau, m3) print(mat.shape) # (rows, 3) import matplotlib.pyplot as plt fig plt.figure(figsize(8, 6)) ax fig.add_subplot(111, projection3d) ax.plot(mat[:, 0], mat[:, 1], mat[:, 2], lw0.5, colorsteelblue) # 三個(gè)軸按實(shí)際數(shù)據(jù)范圍等比防止圖形被壓扁 ax.set_box_aspect((np.ptp(mat[:, 0]), np.ptp(mat[:, 1]), np.ptp(mat[:, 2]))) ax.set_xlabel(x(t)) ax.set_ylabel(x(tτ)) ax.set_zlabel(x(t2τ)) ax.view_init(elev20, azim45) plt.show()邏輯說(shuō)明psr_reconstruct 返回 rows×3 矩陣第 0 列是原序列第 1 列滯后 12 個(gè)采樣點(diǎn)第 2 列滯后 24 個(gè)。等價(jià)于從第 0 個(gè)原始點(diǎn)開(kāi)始以 τ 為步長(zhǎng)取三個(gè)元素構(gòu)成第一個(gè)三維向量然后逐點(diǎn)滑動(dòng)。畫圖用 plot 而不是 scatter幾千個(gè)點(diǎn)只有在連線模式下才能看到連續(xù)的折疊結(jié)構(gòu)線寬 0.5 避免蝶翼兩側(cè)互相糊成一片。set_box_aspect 是三維圖不被壓扁的關(guān)鍵很多流傳的源碼包里沒(méi)有這一句蝴蝶會(huì)被硬拉成飛餅。view_init 固定視角后面做參數(shù)對(duì)比時(shí)才不會(huì)換一個(gè)角度就換了一張圖。mat 行數(shù)超過(guò)兩萬(wàn)時(shí)先隔點(diǎn)抽樣再畫mat[::2] 丟一半點(diǎn)速度翻倍且視覺(jué)幾乎不變。這也是源碼包里經(jīng)常出現(xiàn)的處理不是偷數(shù)據(jù)是控制渲染量。3.3 把 τ 的自動(dòng)估計(jì)接進(jìn)主流程tau_corr autocorr_tau(x, stop1.0 / np.e) tau_mi mi_first_min(x, tau_max80, bins32) print(f自相關(guān)法 tau{tau_corr}, 互信息法 tau{tau_mi}) tau tau_mi if tau_mi is not None else tau_corr mat psr_reconstruct(x, tautau, m3) fig.suptitle(fLorenz, tau{tau}, m3, dt0.02)參數(shù)說(shuō)明自相關(guān)和互信息結(jié)果不一致時(shí)我一般先看一眼互信息曲線確認(rèn)第一極小點(diǎn)旁邊沒(méi)有毛刺再?zèng)Q定是否改用 tau_corr。自動(dòng)估計(jì)的 τ 只配當(dāng)起點(diǎn)不配當(dāng)標(biāo)準(zhǔn)答案——用下一章的參數(shù)掃描驗(yàn)證過(guò)才算數(shù)。4. 相空間重構(gòu)常見(jiàn)問(wèn)題排查五個(gè)翻車現(xiàn)場(chǎng)的現(xiàn)象、原因與對(duì)策相空間重構(gòu)的坑都很隱蔽因?yàn)槌绦虿粫?huì)報(bào)錯(cuò)“τ 選錯(cuò)了”。下面五條按出現(xiàn)頻率排序每一條都值得在自己數(shù)據(jù)上對(duì)照一遍。4.1 現(xiàn)象重構(gòu)軌跡全部擠在空間對(duì)角線附近三維圖是一條細(xì)長(zhǎng)的對(duì)角線或者緊緊貼在一個(gè)平面上看不到蝴蝶的折疊。這是最典型的翻車現(xiàn)場(chǎng)。原因有二τ 太小三個(gè)坐標(biāo)分量數(shù)值幾乎相等或信號(hào)未去均值、帶趨勢(shì)趨勢(shì)項(xiàng)把軌跡拉成一條斜線。經(jīng)驗(yàn)法則凡是吸引子看起來(lái)像個(gè)棒子先懷疑 τ再懷疑預(yù)處理。解決先做預(yù)處理再去調(diào) τ。from scipy.signal import detrend x_clean detrend(x - x.mean()) tau_new mi_first_min(x_clean, tau_max80, bins32) mat psr_reconstruct(x_clean, tau_new, m3)邏輯說(shuō)明detrend 默認(rèn)去掉線性趨勢(shì)去均值消掉直流分量。對(duì)緩慢漂移的實(shí)測(cè)信號(hào)這兩步有時(shí)比調(diào) τ 更關(guān)鍵。處理完再跑互信息法τ 往往會(huì)變大一點(diǎn)軌跡也會(huì)從對(duì)角線上“松開(kāi)”。4.2 現(xiàn)象改變視角后吸引子結(jié)構(gòu)完全變樣同一份數(shù)據(jù)elev20 時(shí)看是蝴蝶elev70 時(shí)看成一團(tuán)亂線兩個(gè)人截圖對(duì)比得出的結(jié)論完全相反。原因三維圖本質(zhì)是二維投影視角和坐標(biāo)縮放都會(huì)扭曲視覺(jué)結(jié)構(gòu)。尤其缺了 set_box_aspect 時(shí)三個(gè)軸按各自范圍獨(dú)立拉伸真實(shí)幾何比例被破壞。解決固定視角加等比盒子。檢查繪圖代碼里有沒(méi)有 set_box_aspect 和 view_init 兩行沒(méi)有就補(bǔ)上。所有參數(shù)對(duì)比統(tǒng)一用同一視角保存圖片時(shí)把視角參數(shù)寫進(jìn)文件名否則截圖無(wú)法追溯。這是血淚經(jīng)驗(yàn)看吸引子形狀必須先固定視角否則等于看圖猜謎。4.3 現(xiàn)象數(shù)據(jù)截?cái)嗪笪咏Y(jié)構(gòu)劇變用前一半數(shù)據(jù)畫圖是一個(gè)環(huán)用后一半畫是另一個(gè)環(huán)掐頭去尾再看形狀大變。原因數(shù)據(jù)里混入了瞬態(tài)段或者系統(tǒng)狀態(tài)本身發(fā)生了遷移。Lorenz 測(cè)試信號(hào)里常見(jiàn)的是初值飛線真實(shí)傳感器數(shù)據(jù)里常見(jiàn)的是緩慢漂移造成的狀態(tài)切換。解決先定位瞬態(tài)段丟掉再用滑動(dòng)窗口截取穩(wěn)態(tài)段。粗略判斷穩(wěn)態(tài)的辦法是計(jì)算每 200 點(diǎn)窗口的質(zhì)心質(zhì)心在三維空間里的偏移超過(guò)坐標(biāo)范圍的 10%就得重新選段。def check_stationary(mat, win200, ratio0.1): center mat.mean(axis0) spans np.ptp(mat, axis0) for start in range(0, len(mat) - win, win): seg_center mat[start:startwin].mean(axis0) if np.any(np.abs(seg_center - center) / spans ratio): return False return True邏輯說(shuō)明質(zhì)心漂移是吸引子結(jié)構(gòu)不穩(wěn)的直接信號(hào)。返回 False 時(shí)別急著調(diào) τ先換數(shù)據(jù)段。這個(gè)函數(shù)對(duì)真實(shí)信號(hào)尤其有用它能直接指出哪一段不屬于同一個(gè)動(dòng)力學(xué)狀態(tài)。4.4 現(xiàn)象τ 選太大軌跡變成稀疏點(diǎn)云三維圖是一堆懸浮的散點(diǎn)看不出連續(xù)軌道像噪聲而非吸引子。原因互信息法自動(dòng)選 τ 時(shí)取錯(cuò)了極小點(diǎn)常見(jiàn)的是第一極小不明顯、代碼誤取第二極小或者信號(hào)周期性太強(qiáng)自相關(guān)法的 1/e 準(zhǔn)則直接失效。解決把互信息曲線畫出來(lái)人工確認(rèn)第一個(gè)極小點(diǎn)。import matplotlib.pyplot as plt taus np.arange(1, 80) mis [mutual_information(x, t, bins32) for t in taus] plt.plot(taus, mis) for i in range(1, len(mis) - 1): if mis[i] mis[i - 1] and mis[i] mis[i 1]: print(局部極小 tau , i 1) plt.show()邏輯說(shuō)明互信息函數(shù)單個(gè) τ 的復(fù)雜度是 O(bins2)80 個(gè) τ 跑下來(lái)也就幾十毫秒放心循環(huán)?？吹角€上低于均值的第一處凹陷那個(gè)位置才是合理 τ不是整條曲線的最低點(diǎn)。如果曲線第一個(gè)極小出現(xiàn)在 tau1說(shuō)明數(shù)據(jù)可能本身采樣過(guò)密或周期性過(guò)強(qiáng)先降采樣再重構(gòu)。4.5 現(xiàn)象兩次運(yùn)行結(jié)果的坐標(biāo)范圍不一致無(wú)法對(duì)比昨天畫的吸引子范圍是 [-20, 20]今天變成 [-15, 15]形狀看著也不一樣但代碼一行沒(méi)改。原因數(shù)據(jù)段起點(diǎn)變了、去趨勢(shì)的位置變了、τ 變了圖上卻看不出參數(shù)差異。這不是算法錯(cuò)誤是復(fù)現(xiàn)管理問(wèn)題。解決每次重構(gòu)輸出時(shí)記錄數(shù)據(jù)段起止索引、τ、m、坐標(biāo)范圍。常見(jiàn)做法是存一個(gè) JSON或者直接編進(jìn)文件名。具體模板放在最后一章這里先記住結(jié)論沒(méi)有參數(shù)快照的重構(gòu)結(jié)果等于沒(méi)有刻度尺的圖紙。5. 三維相空間重構(gòu)的下游定量分析從看圖到算數(shù)三維圖只能讓你“看著像”要說(shuō)服別人、要落到項(xiàng)目里得把“像蝴蝶”變成“D2≈2.05”這種可復(fù)現(xiàn)的數(shù)值。這章講最常用的兩步。5.1 關(guān)聯(lián)維數(shù)G-P 算法把吸引子形狀變成一條飽和曲線from scipy.spatial.distance import pdist def correlation_integral(mat, r): n mat.shape[0] if n 8000: idx np.random.choice(n, 8000, replaceFalse) mat mat[idx] n 8000 dists pdist(mat, metriceuclidean) pairs np.sum(dists r) return 2.0 * pairs / (n * (n - 1))邏輯說(shuō)明pdist 的復(fù)雜度是 O(n2)幾萬(wàn)點(diǎn)會(huì)直接吃爆內(nèi)存所以超過(guò) 8000 行先隨機(jī)抽樣。這里抽的是重構(gòu)軌跡的行也就是相空間里的點(diǎn)不影響幾何結(jié)構(gòu)只降低精度。r 的掃描用對(duì)數(shù)等分rs np.geomspace(0.01, 50, 40) mat3 psr_reconstruct(x, tau, m3) cs np.array([correlation_integral(mat3, r) for r in rs]) # 無(wú)標(biāo)度區(qū)經(jīng)驗(yàn)范圍C(r) 在 0.01 到 0.5 之間 mask (cs 0.01) (cs 0.5) d2 np.polyfit(np.log(rs[mask]), np.log(cs[mask]), 1)[0] print(fD2 ≈ {d2:.3f})參數(shù)說(shuō)明mask 選的是 C(r) 在 0.01 到 0.5 之間的點(diǎn)太小的 r 區(qū)域是離散點(diǎn)噪聲太大則進(jìn)入飽和段。Lorenz 的 D2 文獻(xiàn)值約 2.05算出來(lái)在 1.9 到 2.2 之間都算正常。如果差得遠(yuǎn)不要懷疑算法回去查 τ 和數(shù)據(jù)長(zhǎng)度——這是祖?zhèn)鞯恼{(diào)參順序。5.2 用重構(gòu)軌跡做狀態(tài)識(shí)別的兩個(gè)特征工程落地時(shí)很多人不關(guān)心 D2只想要一個(gè)能區(qū)分“正常”和“異?！钡奶卣?。三維重構(gòu)軌跡可以抽出幾個(gè)比時(shí)域統(tǒng)計(jì)量更敏感的特征。def psr_features(mat): cov np.cov(mat.T) eig np.linalg.eigvalsh(cov) var_ratio np.max(eig) / np.sum(eig) # 主方向方差占比 seg np.diff(mat, axis0) arc_len np.sum(np.linalg.norm(seg, axis1)) # 軌跡總弧長(zhǎng) volume np.prod(np.ptp(mat, axis0)) # 軌跡占據(jù)的空間體積 return var_ratio, arc_len, volume邏輯說(shuō)明var_ratio 反映軌跡在三維空間里鋪得廣不廣結(jié)構(gòu)越扁此值越高arc_len 是軌道在吸引子上繞的總長(zhǎng)度數(shù)據(jù)段相同長(zhǎng)度時(shí)反映繞圈密度volume 是三個(gè)軸范圍的乘積粗估吸引子占據(jù)空間大小。這三個(gè)量對(duì)狀態(tài)切換比均值方差敏感得多。常見(jiàn)做法正常工況取一段數(shù)據(jù)算一組特征異常工況取另一段算一組喂給閾值判斷或 SVM。但要注意邊界特征對(duì)數(shù)據(jù)長(zhǎng)度和預(yù)處理極其敏感對(duì)比時(shí)必須用相同的數(shù)據(jù)段長(zhǎng)度和相同的 τ。比如旋轉(zhuǎn)機(jī)械的振動(dòng)信號(hào)轉(zhuǎn)速一變特征整體漂移得先按轉(zhuǎn)速分段再對(duì)每段單獨(dú)重構(gòu)。5.3 參數(shù)掃描τ 從 1 到 30m 從 3 到 6哪個(gè)組合最穩(wěn)看單張三維圖選 τ 還是容易犯主觀更可靠的辦法是跑參數(shù)掃描看 D2 對(duì)參數(shù)的穩(wěn)定性。results [] for m in [3, 4, 5, 6]: for tau in range(1, 31): mat_t psr_reconstruct(x, tau, mm) rs_t np.geomspace(0.01, 50, 30) cs_t np.array([correlation_integral(mat_t, r) for r in rs_t]) mask_t (cs_t 0.01) (cs_t 0.5) if mask_t.sum() 3: continue d2_t np.polyfit(np.log(rs_t[mask_t]), np.log(cs_t[mask_t]), 1)[0] results.append((m, tau, d2_t))參數(shù)說(shuō)明這組循環(huán)是 4×30120 次 G-P 計(jì)算每次抽樣 8000 點(diǎn)普通筆記本幾分鐘內(nèi)能跑完。選出 D2 隨 m 飽和、且對(duì) τ 變化不敏感的區(qū)域那個(gè) τ 就是穩(wěn)定工作點(diǎn)?！皩?duì) τ 不敏感”本身就是重要信號(hào)——如果 D2 隨 τ 劇烈抖動(dòng)說(shuō)明數(shù)據(jù)長(zhǎng)度不足或系統(tǒng)根本不是單個(gè)吸引子繼續(xù)調(diào)參數(shù)沒(méi)有意義。注意無(wú)標(biāo)度區(qū)的 mask 范圍0.01~0.5只在數(shù)據(jù)量足夠時(shí)有效。數(shù)據(jù)少于 1000 點(diǎn)時(shí)不要強(qiáng)行算 D2結(jié)果沒(méi)有統(tǒng)計(jì)意義。6. 給重構(gòu)結(jié)果留個(gè)狀態(tài)快照文件名就是后悔藥6.1 參數(shù)快照模板與自解釋命名寫完圖或算出 D2 后第一件事是把參數(shù)固化下來(lái)。τ12、m3 這個(gè)組合到底對(duì)應(yīng)哪段數(shù)據(jù)、采樣間隔多少、視角多少度沒(méi)有這些三維圖只是張無(wú)法復(fù)現(xiàn)的插圖。meta { source: lorenz_x, start_idx: 2000, end_idx: 6000, dt: 0.02, tau: 12, m: 3, elev: 20, azim: 45, range: [float(mat.min()), float(mat.max())], } import json with open(recon_meta.json, w) as f: json.dump(meta, f, indent2)參數(shù)說(shuō)明range 記錄三個(gè)軸合并后的最小最大值再次繪圖時(shí)用它統(tǒng)一坐標(biāo)范圍。文件名用“tau12_m3_i2000_6000.png”這種自解釋命名比“重構(gòu)結(jié)果.png”強(qiáng)得多。JSON 里再存一份完整參數(shù)圖丟了還能重建。6.2 換數(shù)據(jù)前的內(nèi)置校驗(yàn)我被這類問(wèn)題坑過(guò)不止一次同一份振動(dòng)數(shù)據(jù)上午下午各跑一遍畫出的圖一個(gè)寬一個(gè)扁最后發(fā)現(xiàn)只是一個(gè) τ 用 8、一個(gè)用 10還沒(méi)人記得誰(shuí)用了哪個(gè)。從那以后所有重構(gòu)實(shí)驗(yàn)一律帶參數(shù)快照。一個(gè)實(shí)用的驗(yàn)證習(xí)慣把代碼換到陌生數(shù)據(jù)上之前先在 Lorenz 上復(fù)現(xiàn) D2≈2.05確認(rèn)整個(gè)代碼通道沒(méi)問(wèn)題再碰真實(shí)數(shù)據(jù)。真實(shí)數(shù)據(jù)算出的 D2 落在 1.1 到 2.9 之間通常說(shuō)明有確定性結(jié)構(gòu)接近整數(shù)或半整數(shù)更有說(shuō)服力如果 D2 大于 4 或找不到無(wú)標(biāo)度區(qū)先懷疑數(shù)據(jù)而不是算法。真正常規(guī)、能反復(fù)用、能對(duì)比的相空間重構(gòu)流程一定長(zhǎng)著“參數(shù)看得見(jiàn)、視角固定、坐標(biāo)等比”的樣子。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取