【脑电6】

脑电 6

目录

• 1. 为什么需要非线性分析
• 2. 熵——信号不可预测性的度量
• 3. 复杂度——信号中模式丰富度的量化
• 4. 分形与自相似性
• 5. 去趋势波动分析——长程相关性
• 6. 非线性特征工具箱

1. 为什么需要非线性分析

1.1 脑电不是线性系统

线性系统的典型特征：输入和输出成比例，不同频率成分互不干扰。但大脑是典型的非线性系统：

• 两个刺激同时呈现 ≠ 单独呈现的EEG响应之和
• α波和γ波之间存在跨频率耦合（γ功率随α相位变化）
• 癫痫从正常→发作的转变是突然的相变（非线性动力学中的分岔现象）
• 同一被试在"相同"状态下的EEG波形从不完全重复

importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportsignal# 直观演示：线性系统 vs 非线性系统t=np.linspace(0,4,500)x=np.sin(2*np.pi*2*t)+0.5*np.sin(2*np.pi*8*t)# 线性系统：y = 2*x（简单放大）y_linear=2*x# 非线性系统：y = x³（输出不是输入的简单缩放）y_nonlinear=x**3fig,axes=plt.subplots(2,2,figsize=(14,7))axes[0,0].plot(t,x,'b',lw=1);axes[0,0].set_title('输入信号 x(t)')axes[0,1].plot(t,y_linear,'g',lw=1);axes[0,1].set_title('线性系统输出 y=2x\n(形状不变)')axes[1,0].plot(t,y_nonlinear,'r',lw=1);axes[1,0].set_title('非线性系统输出 y=x³\n(波形变形、高频增加)')# 输入-输出散点图axes[1,1].scatter(x,y_linear,s=2,alpha=0.3,label='线性(y=2x)')axes[1,1].scatter(x,y_nonlinear,s=2,alpha=0.3,label='非线性(y=x³)')axes[1,1].set_xlabel('输入');axes[1,1].set_ylabel('输出')axes[1,1].legend();axes[1,1].set_title('输入-输出关系')plt.tight_layout();plt.show()print("线性：输入×2 = 输出×2，比例恒定")print("非线性：输入和输出的关系不是直线——这就是非线性分析的动机")

1.2 非线性指标能捕捉的东西

2. 熵——信号不可预测性的度量

2.1 三种熵的对比

defsample_entropy(data,m=2,r_factor=0.2):""" 样本熵——衡量信号中新模式产生的概率 熵高=信号不可预测=混沌/警觉 熵低=信号规律=有序/冥想/麻醉 """N,r=len(data),r_factor*np.std(data)defcount_matches(template_len):count=0templates=np.array([data[i:i+template_len]foriinrange(N-template_len)])foriinrange(len(templates)):dists=np.max(np.abs(templates-templates[i]),axis=1)count+=np.sum(dists<r)-1# 减1排除自身匹配returncount A,B=count_matches(m+1),count_matches(m)return-np.log(A/(B+1e-10))ifB>0elsefloat('inf')defpermutation_entropy(data,m=3,delay=1):""" 排列熵——只看数值大小排序模式 例：窗口[1,5,3]→排序[0,2,1]→"中-小-大"模式 对所有窗口统计模式分布→算熵 """fromitertoolsimportpermutations N=len(data)patterns=list(permutations(range(m)))pattern_count={p:0forpinpatterns}total=0foriinrange(N-(m-1)*delay):window=data[i:i+m*delay:delay]# 找到排序模式sorted_idx=tuple(np.argsort(window))ifsorted_idxinpattern_count:pattern_count[sorted_idx]+=1total+=1# 计算熵pe=0forcountinpattern_count.values():ifcount>0:p=count/total pe-=p*np.log(p)returnpe/np.log(len(patterns))# 归一化到[0,1]# 对比演示np.random.seed(42)noise=np.random.randn(1000)# 白噪声→高熵sine=np.sin(2*np.pi*np.arange(1000)/256*10)# 正弦→低熵mixed=noise*0.6+sine*0.4# 混合→中等print("信号类型对比:")forname,sigin[('白噪声(高熵)',noise),('正弦波(低熵)',sine),('混合(中熵)',mixed)]:se=sample_entropy(sig)pe=permutation_entropy(sig[:500])print(f"{name}: SampleEn={se:.3f}, PermEn={pe:.3f}")

3. 复杂度——信号中模式丰富度的量化

3.1 Lempel-Ziv复杂度

LZC衡量的不是"信号有多乱"，而是"信号中有多少种不同的子串模式"——模式越多=越复杂。

deflempel_ziv_complexity(data):"""LZ复杂度——计算信号二值化后不同子串模式数量"""# 二值化：大于中位数=1，小于=0binary=(data>np.median(data)).astype(int)n=len(binary)i,C,L=0,1,1# C=复杂度计数器, L=当前窗口长度whilei+L<=n:pattern=''.join(map(str,binary[i:i+L]))history=''.join(map(str,binary[:i+L-1]))ifpatterninhistory:L+=1# 模式出现过→扩大窗口else:C+=1# 新模式→复杂度+1i+=L L=1returnC/n# 归一化# 直观演示：规则序列 vs 随机序列regular=np.tile([1,0,1,0],50)# 高度规则→低复杂度random_seq=np.random.randint(0,2,200)# 随机→高复杂度print(f"规则序列(1010重复) LZC:{lempel_ziv_complexity(regular):.3f}")print(f"随机序列 LZC:{lempel_ziv_complexity(random_seq):.3f}")print("→ 冥想深度↑ → 大脑活动模式趋于简单 → LZC↓")

3.2 LZC在脑电中的应用

• 麻醉深度监测：LZC随麻醉加深而持续降低——比频谱指标更可靠
• 冥想状态评估：深度冥想时LZC显著低于静息状态
• 癫痫预测：发作前几分钟LZC异常降低——可用于预警
• 意识障碍评估：植物状态患者的LZC显著低于最小意识状态

4. 分形与自相似性

4.1 脑电是分形信号

分形的核心特征：自相似性——放大看局部，和整体有相似的统计特性。EEG信号在毫秒到分钟的多个时间尺度上都呈现1/f频谱（功率∝1/f^α），这正是分形信号的特征。

defhiguchi_fractal_dimension(data,kmax=10):""" Higuchi分形维度——衡量信号的"曲折程度" D范围：1(光滑规则曲线) ~ 2(极其曲折，填满平面) 冥想加深→D降低（信号更光滑规则） 麻醉加深→D降低 癫痫发作→D降低 """N=len(data)Lk=np.zeros(kmax)forkinrange(1,kmax+1):Lmk=np.zeros(k)forminrange(k):idxs=np.arange(m,N-1,k)iflen(idxs)>1:Lmk[m]=np.sum(np.abs(np.diff(data[idxs])))Lmk[m]*=(N-1)/((len(idxs)-1)*k)valid=Lmk[Lmk>0]Lk[k-1]=np.mean(valid)iflen(valid)>0else0# 最小二乘拟合: log(Lk) ~ log(1/k)x=np.log(1.0/np.arange(1,kmax+1))y=np.log(Lk+1e-10)D=-np.polyfit(x,y,1)[0]returnDdefpetrosian_fractal_dimension(data):"""Petrosian分形维度——基于零交叉点计数，计算极快"""# 零交叉：信号穿过均值的次数zero_crossings=np.sum(np.diff(data>np.mean(data))!=0)N=len(data)returnnp.log10(N)/(np.log10(N)+np.log10(N/(N+0.4*zero_crossings)))

5. 去趋势波动分析——长程相关性

5.1 DFA的原理

DFA（Detrended Fluctuation Analysis）衡量信号在不同时间尺度上的波动特征，核心输出是一个指数α：

• α ≈ 0.5：白噪声——各时间尺度随机，无长程相关
• α ≈ 1.0：1/f噪声（粉红噪声）——健康大脑的典型特征
• α > 1.0：存在长程正相关——前一刻的趋势延续到下一刻
• α < 0.5：反相关——涨落后倾向于被跌落跟随

defdfa(data,scales=None):""" 去趋势波动分析(DFA) 返回缩放指数α——衡量长程相关性强度 """ifscalesisNone:scales=np.unique(np.logspace(1,np.log10(len(data)//4),10).astype(int))# 1. 累积和（积分）y=np.cumsum(data-np.mean(data))# 2. 对每个尺度计算波动fluctuations=np.zeros(len(scales))foridx,scaleinenumerate(scales):n_segments=len(data)//scale rms_total=0foriinrange(n_segments):seg=y[i*scale:(i+1)*scale]# 用多项式拟合去趋势（线性趋势=DFA1）x=np.arange(len(seg))trend=np.polyval(np.polyfit(x,seg,1),x)detrended=seg-trend rms_total+=np.mean(detrended**2)fluctuations[idx]=np.sqrt(rms_total/n_segments)# 3. 拟合log(F) vs log(scale) → 斜率=αvalid=~np.isnan(fluctuations)&(fluctuations>0)alpha=np.polyfit(np.log10(scales[valid]),np.log10(fluctuations[valid]),1)[0]returnalpha,scales,fluctuations# 演示：不同α的对比np.random.seed(42)white_noise=np.random.randn(2000)# α≈0.5# 生成1/f噪声（简化的Pink噪声）pink_noise=np.cumsum(np.random.randn(2000))# α≈1.5（布朗噪声）forname,sigin[('白噪声(α≈0.5)',white_noise),('布朗噪声(α≈1.5)',pink_noise)]:alpha,_,_=dfa(sig)print(f"{name}: DFA α ={alpha:.3f}")print("健康静息EEG的DFA α通常在0.7-1.0之间")print("深度麻醉/昏迷时α可能下降到0.5左右")

5.2 DFA在临床和研究中的应用

6.非线性特征工具箱

classNonlinearEEGAnalyzer:"""非线性EEG特征提取工具箱"""def__init__(self,fs=256):self.fs=fsdefextract_all(self,eeg_segment):"""一键提取所有非线性特征"""features={}# 熵类features['sample_entropy']=sample_entropy(eeg_segment)features['perm_entropy']=permutation_entropy(eeg_segment,m=3)# 复杂度features['lz_complexity']=lempel_ziv_complexity(eeg_segment)# 分形features['higuchi_fd']=higuchi_fractal_dimension(eeg_segment)features['petrosian_fd']=petrosian_fractal_dimension(eeg_segment)# 长程相关alpha,_,_=dfa(eeg_segment)features['dfa_alpha']=alphareturnfeaturesdefcompare_states(self,eeg_segments,state_names):"""对比不同状态的非线性特征"""results=[]forseg,nameinzip(eeg_segments,state_names):results.append(self.extract_all(seg))# 可视化feat_names=list(results[0].keys())n_feat=len(feat_names)fig,axes=plt.subplots(1,n_feat,figsize=(3*n_feat,4))fori,(feat,ax)inenumerate(zip(feat_names,axes)):values=[r[feat]forrinresults]ax.bar(range(len(state_names)),values,color=['steelblue','salmon','orange','green'][:len(values)],edgecolor='black')ax.set_xticks(range(len(state_names)))ax.set_xticklabels(state_names,fontsize=7)ax.set_title(feat,fontsize=9);ax.grid(True,alpha=0.3,axis='y')plt.suptitle('不同脑状态下非线性特征对比',fontsize=13)plt.tight_layout();plt.show()returnresults# 演示：四种模拟状态analyzer=NonlinearEEGAnalyzer(fs=256)np.random.seed(42)t=np.arange(0,4,1/256)# 模拟四种状态rest=20*np.sin(2*np.pi*10*t)+np.random.randn(len(t))*3meditation=25*np.sin(2*np.pi*10*t)+np.random.randn(len(t))*1.5# 更干净规则drowsy=10*np.sin(2*np.pi*5*t)+np.random.randn(len(t))*5# 更多慢波噪声noise_sig=np.random.randn(len(t))*8# 纯噪声results=analyzer.compare_states([rest,meditation,drowsy,noise_sig],['静息','深度冥想','困倦','随机噪声'])print("\n各状态特征汇总:")fori,(name,r)inenumerate(zip(['静息','深度冥想','困倦','随机噪声'],results)):print(f"\n{name}:")fork,vinr.items():print(f"{k}:{v:.4f}")ifi==1:print(" → 深度冥想：熵最低、复杂度最低、分形维最低（最有序）")