丘腦底核神經活動與帕金森病癥狀密切相關,其局部場電位提供了豐富的神經波動成分信息。在震顫狀態下,局部場電位信號中可見震顫頻率及其倍頻譜峰,這種多倍頻現象中各成分間的相互關系很可能與震顫發生的神經機制有關。為確定各神經活動成分間的相互耦合關系,對9名帕金森病患者局部場電位進行跨頻率相關分析,結果表明在震顫頻率與其倍頻成分之間存在明顯相關,互相關系數在0.18~0.50范圍內。Granger因果分析表明震顫基頻對2倍頻成分、震顫基頻和2倍頻成分對3倍頻的發生有貢獻作用。在雙譜分析中發現震顫基頻和2倍頻、3倍頻成分間存在二次相位耦合作用,說明高頻段倍頻成分的產生與低頻段倍頻成分間的非線性耦合有關。本研究證實了震顫相關的基底神經節活動的多個波動成分并不獨立,互相間存在線性和非線性交互作用關系,進一步揭示了震顫發生的復雜神經機制。
引用本文： 王宗寶, 黃永志, 張新靜, 耿馨佚, 陳驍, 王守巖. 帕金森病患者局部場電位信號多頻耦合特征分析. 生物醫學工程學雜志, 2015, 32(4): 874-880. doi: 10.7507/1001-5515.20150156 復制
引言
帕金森病患者丘腦底核的異常神經活動與疾病癥狀密切相關[1]，其中過度同步化的beta頻段(12~30 Hz)活動是導致帕金森病患者行動遲緩、肌肉強直等癥狀的重要因素[2]，而電刺激抑制beta活動可以達到緩解這些癥狀的效果[3]。但丘腦底核神經活動在震顫癥狀發生中的作用仍未明確。僅在部分患者可見丘腦底核震顫頻率處的神經活動隨震顫發生顯著增加[4-5]，可能是震顫發作與多個神經活動成分的整合而非其絕對幅度有關[5]。
帕金森病患者丘腦底核局部場電位的0~20 Hz頻段，有明顯的類似諧波的頻譜峰值[6-7]，例如4.5、9.0、13.5和18.0 Hz處譜峰，這些以震顫頻率為基頻的倍頻成分間的相互關系對解釋震顫發生機制有重要意義。以震顫頻率和其2倍頻成分為例，對震顫發作時的場電位和肌電信號的相干譜分析顯示，兩種信號在這兩個頻率處均有顯著的相干性[7]，但無法確定2倍震顫頻率處的成分為生理信號，還是震顫頻率成分的二次諧波。如果震顫頻率成分和其2倍頻成分互相獨立，則表明神經環路中存在兩個獨立的振子，震顫的發生可能是兩個振子耦合的結果，這有可能是震顫不穩定的原因。
多個倍頻成分存在有兩種假設：一種假設是成分間互相獨立，均源于自發的神經活動，由神經活動的內在機制產生，每個成分對應了大腦活動的一種特定模式，不同成分的特征可用于對大腦活動狀態進行判斷，例如12~30 Hz頻段能量異常瞬態增加是帕金森病發病的一個重要病理特征[8]；另一種假設是，局部場電位不同頻率成分之間可能存在著線性或非線性耦合關系，而多倍頻成分現象是耦合作用的外在信號特征。
局部場電位成分間的線性耦合關系可基于信號幅度的相關分析和因果分析進行研究。首先利用功率譜分析確定局部場電位信號頻譜譜峰的分布，以及利用短時傅里葉變換(short time Fourier transform,STFT)方法可觀察不同頻率成分的能量隨時間的動態變化[9]。其次可利用相關分析對相應頻率成分的幅度包絡進行線性相關程度度量。進而可利用Granger因果分析度量一個頻率成分的幅度對另一個成分幅度的改變是否存在因果性，考察基頻與倍頻成分幅度間的雙向線性因果關系[10]。
局部場電位各成分間的非線性耦合關系可用高階譜方法分析觀察。神經活動通常具有結構復雜和信息豐富的特點，腦深部局部場電位不同成分可能通過某些非線性過程對其他頻率成分產生影響，雙譜分析可檢測非平穩信號內部的二次相位耦合機制[11-12]，這種非線性作用是指兩個頻率成分在滿足頻率耦合和相位耦合條件時可激發產生第三個頻率成分。使用高階統計量中的雙譜可以檢測二次相位耦合是否存在以及度量耦合強度，因而可以確定局部場電位各成分之間的二次非線性關系。
本文為研究丘腦底核局部場電位多倍頻成分間的相互作用關系，采用功率譜分析和時頻分析方法觀察帕金森病患者局部場電位各頻段的能量特征和各頻率成分隨時間的動態變化；采用相關分析、Granger因果分析方法，度量多倍頻成分的幅度包絡間的線性相關程度和雙向線性因果關系；采用雙譜分析方法檢測局部場電位內部二次相位耦合作用。
1 方法
1.1 局部場電位記錄與預處理
共有9名帕金森病患者參與本實驗，在英國牛津大學John Radcliffe Hospital完成丘腦底核深部腦刺激電極植入手術，植入電極為美敦力DBS 3389，電極直徑為1.27 mm，各電極觸點長度為1.5 mm，4個觸點之間的相鄰間隔為0.5 mm，電極植入后通過核磁共振成像對位置進行確認。
電極植入4~6 d后，從外置的電極導線采用雙極方法(相鄰觸點0-1,1-2,2-3)記錄丘腦底核深部局部場電位。從震顫相關的前臂伸肌和屈肌肌群同步記錄表面肌電圖。信號采樣頻率均為2 000 Hz，信號經由CED1902隔離放大器放大(Cambridge Eelectronic Design，UK)，局部場電位放大10 000倍，肌電信號放大1 000倍，帶通濾波參數為0.5~500 Hz，并使用CED1401 II和Spike2軟件(Cambridge Eelectronic Design，UK)對信號進行模數轉換和采集。實驗獲得當地道德倫理委員會批準許可。數據已得到英國牛津大學功能神經外科研究組授權使用。
預處理對信號進行90 Hz抗混疊低通濾波(切比雪夫I型濾波器)，降采樣2 000 Hz信號至250 Hz，采用自適應迭代濾波除去50 Hz工頻干擾，最后進行3 Hz高通濾波去除低頻成分(切比雪夫I型濾波器)。
1.2 時頻分析
應用STFT方法觀察局部場電位不同成分的能量密度隨時間的變化情況。離散隨機信號x(i),i=1,2,3…的STFT定義為
| ${\rm{STFT}}\left( {k,m} \right) = \sum\limits_{i = 1}^L {x\left( i \right){h^ * }\left( {i - k} \right)W_L^{mi}} ,$ |
其中窗函數h(k)時刻k為中心，每一個時刻k都可以得到一個不同的頻譜，則時頻譜密度P(k,ω)為
| $P\left( {k,\omega } \right) = {\left| {{\rm{STFT}}\left( {k,m} \right)} \right|^2}$ |
1.3 跨頻率相關分析
根據不同患者信號的譜峰分布，采用帶通濾波提取震顫相關的頻率成分，通過Hilbert變換提取各成分的幅度包絡，進而用互相關函數定量分析各成分包絡信號之間的相關程度。對于有限長的樣本，互相關函數的估計為
| ${\hat r_{xy}}\left( m \right) = \frac{1}{N}\sum\limits_{n = 0}^{N - 1 - \left| m \right|} {x\left( n \right)y\left( {n + m} \right)} \;\;0 \le \left| m \right| \le N - 1$ |
其中x(i)和y(i)的長度均為N，互相關函數序列長度為2N-1。
1.4 時域Granger因果分析
在Hilbert變換提取不同頻率成分的幅度包絡信息基礎上，Granger因果分析可對兩組時間序列之間的定向線性因果關系進行度量，衡量一個時間序列x(t)是否對另一個時間序列y(t)的發生有貢獻，進而估計兩個信號之間存在的因果影響的程度。
Granger因果分析簡單定義如下，若兩個廣義平穩時間序列x(t)和y(t)之間的線性關系可以用雙變量AR模型來表示，即
| $\begin{array}{l} x\left( t \right) = \sum\limits_{k = 1}^P {{a_{11}}\left( k \right)x\left( {t - k} \right)} + \sum\limits_{k = 1}^P {{a_{12}}\left( k \right)y\left( {t - k} \right)} + {v_1}\left( t \right)\\ y\left( t \right) = \sum\limits_{k = 1}^P {{a_{21}}\left( k \right)x\left( {t - k} \right)} + \sum\limits_{k = 1}^P {{a_{22}}\left( k \right)y\left( {t - k} \right)} + {v_2}\left( t \right), \end{array}$ |
其中p為模型階數，aij(k)為模型參數，v1(t)和v2(t)分別為x(t)和y(t)的預測誤差。當x(t)和y(t)分別僅用自身作預測，v1(t)和v2(t)的方差分別為Σ1和Σ2。當x(t)和y(t)利用自身和對方的前p個值作預測，此時v1(t)和v2(t)的方差分別為Σ11和Σ22。對于y(t)引起x(t)的因果關系時域上可以定義為
| ${F_{y \to x}} = \ln \frac{{\sum 1 }}{{\sum 1 1}}\left( {{f_{y \to x}} \ge 0} \right)$ |
若Σ11<Σ1，則說明加入y(t)對預測x(t)有幫助，預測誤差減少，認為y(t)對x(t)的發生有貢獻，反之亦然。
1.5 雙譜分析
高階譜分析中的雙譜分析定義為隨機時間信號三階累積量的二維傅里葉變換：
| $\begin{array}{l} R\left( {{m_1},{m_2}} \right) = E\left[ {x\left( n \right)x\left( {n + {m_1}} \right)x\left( {n + {m_2}} \right)} \right]\\ B\left( {{\omega _1},{\omega _2}} \right) = \sum\limits_{{m_1} = - \infty }^\infty {\sum\limits_{{m_2} = - \infty }^\infty {R\left( {{m_1},{m_2}} \right)} } \exp \left[ { - j\left( {{\omega _1}{m_1} + {\omega _2}{m_2}} \right)} \right], \end{array}$ |
其中R(m1,m2)為隨機信號x(t)的三階累積量，m1和m2為時間延遲量，B(ω1,ω2)為頻率對(ω1,ω2)的雙譜值。
雙譜保留了信號的相位信息，可以檢測出兩個頻率成分與它們的和頻率成分之間有無二次相位耦合作用[13]。二次相位耦合指一個非線性作用過程中兩個頻率成分在同時滿足頻率耦合和相位耦合的條件下，會在它們的和頻率處激發產生能量的現象。由下式中兩個互相獨立的隨機信號y(t)和z(t)為例說明：
| $\begin{array}{l} y\left( t \right) = \cos \left( {2\pi {\omega _1}t + {\varphi _1}} \right) + \cos \left( {2\pi {\omega _2} + {\varphi _2}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\cos \left( {2\pi {\omega _3} + {\varphi _3}} \right)\\ z\left( t \right) = \cos \left( {2\pi {\omega _1}t + {\varphi _1}} \right) + \cos \left( {2\pi {\omega _2} + {\varphi _2}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\cos \left( {2\pi {\omega _3} + {\varphi _4}} \right), \end{array}$ |
y(t)和z(t)均滿足ω1+ω2=ω3,φ1,φ2,φ3是[-π,π]內均勻分布的獨立隨機變量，且φ1+φ2=φ4。信號z(t)中存在二次相位耦合作用，ω1,ω2兩個成分的二次相位耦合對ω3處能量的產生有貢獻作用，而信號y(t)中不存在二次相位耦合作用，ω3處能量產生與ω1,ω2成分無關。
采用間接法估計雙譜，選用Parzen窗，每段數據長2 s，時間延遲為1.2 s，數據重疊1 s，1 024點FFT。在實際應用中，受到估計誤差和信號長度的限制，雙譜域中可能出現僅有頻率耦合沒有相位耦合的偽峰。替代數據方法消除了原信號中原有的非線性過程，僅保留了與原信號相同的平均功率譜密度特征，是鑒別偽峰的有效方法[14]。采用AAFT方法產生替代數據[15]，將原數據產生的雙譜與替代數據產生的雙譜進行比較，若兩者之間有顯著性差異(P<0.05，t-test)，則代表有統計學意義的雙譜譜峰存在。
2 結果
2.1 譜分析與時頻分析
在深部腦刺激電極植入后，在震顫和靜息狀態下同步記錄肌電信號和丘腦底核局部場電位，功率譜分析顯示9名患者丘腦底核場電位信號在震顫狀態下均存在一定程度的似諧波現象，如表 1所示，所有患者均存在震顫頻率(基頻)的譜峰(3~5 Hz)和2倍頻譜峰，5名患者中存在3倍頻譜峰。
表1
9例帕金森病患者震顫狀態丘腦底核局部場電位譜峰分布
Table1.
Distribution of power spectral peaks in tremor LFPs recorded from STNs of 9 Parkinson’s disease patients
以#1的40 s的肌電信號和局部場電位信號，如圖 1(a)、(b)所示為例說明，0~28 s時間段帕金森病患者處于震顫狀態，28~40 s時間段帕金森病患者處于靜息狀態。震顫狀態下3.5、7.0、10.5 Hz處有明顯的功率譜峰，如圖 1(c)所示；靜息狀態沒有這種現象，如見圖 1(d)所示。圖 1(e)為#1局部場電位信號的STFT結果，在震顫狀態時，4個頻率成分均保持較高的能量密度，28 s后從震顫狀態轉為靜息狀態，10.5 Hz和14.0 Hz能量密度急劇下降，而3.5 Hz和7.0 Hz處能量密度則在32 s后出現減弱，在其余患者均可見震顫狀態下各倍頻成分能量密度高于靜息狀態，以及各譜峰幅度隨時間而動態變化的現象，進而應用相關分析和因果分析度量各譜峰幅度動態變化之間的關系。
圖1
#1局部場電位時頻分析
(a)震顫狀態(0~28 s)和靜息狀態(28~40 s)的相關肌電信號；(b)與A同步記錄的震顫-靜息狀態的丘腦底核局部場電位；(c)震顫狀態丘腦底核局部場電位平均功率譜密度；(d)靜息狀態丘腦底核局部場電位平均功率譜密度；(e) 是(b)中局部場電位的時頻譜圖
Figure1. Time-frequency analysis of LFP signals in #1 patientTremor (0~28 s) and resting (28~40 s) related EMGs (a),LFPs (b),power spectrum of tremor (c) and resting (d) states,short-time Fourier transform of LFPs (e)
2.2 局部場電位各成分相關分析
對9例帕金森病患者震顫狀態的場電位信號進行跨頻率相關分析。根據每例患者局部場電位信號譜峰分布位置不同，分別對各例患者信號進行帶通濾波，提取各基頻、2倍頻和3倍頻頻段成分。應用Hilbert變換提取各成分的幅度包絡，計算兩兩間的互相關函數。圖 2為9名患者各頻率成分包絡間的最大互相關系數，倍頻成分的包絡信號互相關系數在0.18~0.50范圍內，說明各倍頻成分的幅度包絡信號間保持一定的相關性。
圖2
震顫狀態下9例帕金森病患者各成分包絡之間的平均最大互相關系數(均值±標準差)
Figure2.
Average of maximum cross-correlation coefficients between different components in LFPs of 9 Parkinson’s disease patients during tremor state (mean±standard deviation)
2.3 局部場電位各成分Granger因果分析
應用帶通濾波和Hilbert變換得到場電位各頻率成分的包絡，進而可利用時域Granger因果分析的方法估計各成分間的因果關系，平均因果分析結果如圖 3所示。9例帕金森病患者信號的實驗結果平均值顯示，震顫基頻對2倍頻(顯著性檢驗P<0.05，paired t-test)、震顫基頻對3倍頻、2倍頻對3倍頻有較強的因果值，說明處于較高頻段的倍頻成分的產生受到處于低頻段的基頻或倍頻成分的影響較大。
圖3
9名帕金森病患者震顫狀態下場電位各成分包絡之間的Granger因果分析(均值±標準差)
Figure3.
Granger causality between different components in LFPs of 9 Parkinson’s disease patients during tremor state (mean±standard deviation)
2.4 局部場電位雙譜分析
對震顫和靜息狀態下的局部場電位進行雙譜分析和統計學檢驗。9例帕金森病患者局部場電位信號在震顫狀態下的二次相位耦合作用主要集中基頻與基頻、2倍頻、3倍頻之間，如圖 4所示，說明在震顫狀態下震顫頻率自身耦合對2倍頻處的能量峰值產生有貢獻作用，震顫頻率和2倍頻之間的二次相位耦合對3倍頻處的能量產生有貢獻作用。這證明了震顫頻率相關的神經波動的存在，以及與其他倍頻成分的耦合作用驅動了更高頻率，如beta頻段成分的產生。靜息狀態下非線性耦合作用在雙譜域上的分布位置則不能解釋倍頻成分之間的驅動關系。
圖4
9例患者震顫狀態(實心圓)和靜息狀態(空心圓)信號雙譜譜峰位置分布
Figure4.
Bispectral peak positions in tremor state (filled circles) and resting state (unfilled circles) of 9 Parkinson’s disease patients
以患者#1為例進一步說明。與替代數據結果比較進行統計學顯著性檢測，患者#1震顫狀態下的場電位雙譜分析結果表明有6個譜峰如圖 5(a)、(c)所示，分別位于雙譜域中基頻-基頻、基頻-2倍頻、基頻-3倍頻的相應位置(P<0.05)，而在靜息狀態下僅在基頻-基頻處有一相對減弱的譜峰(P<0.05)。這說明在震顫狀態下，震顫基頻的自身耦合對2倍頻處的能量峰值產生有貢獻作用，震顫基頻和倍頻之間的二次相位耦合對3倍頻處的能量產生有貢獻作用。
圖5
局部場電位雙譜分析結果
#1震顫狀態丘腦底核局部場電位的雙譜幅度(a)和譜峰位置(c)；#1靜息狀態丘腦底核局部場電位的雙譜幅度(b)和譜峰位置(d)
Figure5. Bispectral analysis of LFP signalsAmplitude (a) and peak positions (c) of bispectra in tremor state of #1; Amplitude (b) and peak positions (d) of bispectra in resting state of #1
3 討論與結論
本文研究了帕金森病患者丘腦底核局部場電位在震顫和靜息兩種帕金森病疾病狀態下神經波動成分之間的作用關系，在患者震顫狀態的局部場電位功率譜中發現了“m×震顫頻率”的似諧波現象(m=1,2,3)。進而應用跨頻率相關分析、時域Granger因果分析和雙譜分析，研究了各倍頻成分之間的線性和非線性作用關系。
時頻分析結果顯示，震顫基頻和2倍頻成分集中了震顫狀態主要的能量分布，基頻、2倍頻成分和3倍頻保持一定的同步性。跨頻率相關分析發現，各倍頻成分的幅度間保持0.18~0.50的相關性，表明各成分間存在相關關系，一定程度上解釋了信號時頻譜圖中能量密度變化的現象。時域Granger因果分析量化了一個成分對產生另一個成分的影響，
實驗結果顯示了低頻處的基頻和2倍頻成分對高頻處的3倍頻成分相比高頻成分對低頻成分有更強的線性因果關系。
雙譜分析檢測了倍頻成分間的二次非線性關系，解釋了震顫基頻自身、震顫基頻與2倍頻頻段倍頻成分因為相位耦合作用激發了2倍頻、3倍頻成分的能量，進而解釋了時頻分析中震顫到靜息狀態轉變過程中各成分能量密度的變化現象。
震顫是帕金森病的典型癥狀，其表現形式為有規律的肢體顫動，藥物治療對震顫的控制作用較為有限。研究人員對震顫的發病原因提出過多種假設，但至今其神經機制仍不清楚。本工作在我們以往的工作基礎上進一步證實震顫的產生與基底神經節神經活動的多個波動成分有關，而且這些成分并非獨立，存在著相互之間的線性和非線性交互作用，并具有隨時間改變的動態性，且與患者個體有關，而這些神經活動是由一個神經環路亦或是來自于多個神經環路匯集于基底神經節的仍待進一步探索。
引言
帕金森病患者丘腦底核的異常神經活動與疾病癥狀密切相關[1]，其中過度同步化的beta頻段(12~30 Hz)活動是導致帕金森病患者行動遲緩、肌肉強直等癥狀的重要因素[2]，而電刺激抑制beta活動可以達到緩解這些癥狀的效果[3]。但丘腦底核神經活動在震顫癥狀發生中的作用仍未明確。僅在部分患者可見丘腦底核震顫頻率處的神經活動隨震顫發生顯著增加[4-5]，可能是震顫發作與多個神經活動成分的整合而非其絕對幅度有關[5]。
帕金森病患者丘腦底核局部場電位的0~20 Hz頻段，有明顯的類似諧波的頻譜峰值[6-7]，例如4.5、9.0、13.5和18.0 Hz處譜峰，這些以震顫頻率為基頻的倍頻成分間的相互關系對解釋震顫發生機制有重要意義。以震顫頻率和其2倍頻成分為例，對震顫發作時的場電位和肌電信號的相干譜分析顯示，兩種信號在這兩個頻率處均有顯著的相干性[7]，但無法確定2倍震顫頻率處的成分為生理信號，還是震顫頻率成分的二次諧波。如果震顫頻率成分和其2倍頻成分互相獨立，則表明神經環路中存在兩個獨立的振子，震顫的發生可能是兩個振子耦合的結果，這有可能是震顫不穩定的原因。
多個倍頻成分存在有兩種假設：一種假設是成分間互相獨立，均源于自發的神經活動，由神經活動的內在機制產生，每個成分對應了大腦活動的一種特定模式，不同成分的特征可用于對大腦活動狀態進行判斷，例如12~30 Hz頻段能量異常瞬態增加是帕金森病發病的一個重要病理特征[8]；另一種假設是，局部場電位不同頻率成分之間可能存在著線性或非線性耦合關系，而多倍頻成分現象是耦合作用的外在信號特征。
局部場電位成分間的線性耦合關系可基于信號幅度的相關分析和因果分析進行研究。首先利用功率譜分析確定局部場電位信號頻譜譜峰的分布，以及利用短時傅里葉變換(short time Fourier transform,STFT)方法可觀察不同頻率成分的能量隨時間的動態變化[9]。其次可利用相關分析對相應頻率成分的幅度包絡進行線性相關程度度量。進而可利用Granger因果分析度量一個頻率成分的幅度對另一個成分幅度的改變是否存在因果性，考察基頻與倍頻成分幅度間的雙向線性因果關系[10]。
局部場電位各成分間的非線性耦合關系可用高階譜方法分析觀察。神經活動通常具有結構復雜和信息豐富的特點，腦深部局部場電位不同成分可能通過某些非線性過程對其他頻率成分產生影響，雙譜分析可檢測非平穩信號內部的二次相位耦合機制[11-12]，這種非線性作用是指兩個頻率成分在滿足頻率耦合和相位耦合條件時可激發產生第三個頻率成分。使用高階統計量中的雙譜可以檢測二次相位耦合是否存在以及度量耦合強度，因而可以確定局部場電位各成分之間的二次非線性關系。
本文為研究丘腦底核局部場電位多倍頻成分間的相互作用關系，采用功率譜分析和時頻分析方法觀察帕金森病患者局部場電位各頻段的能量特征和各頻率成分隨時間的動態變化；采用相關分析、Granger因果分析方法，度量多倍頻成分的幅度包絡間的線性相關程度和雙向線性因果關系；采用雙譜分析方法檢測局部場電位內部二次相位耦合作用。
1 方法
1.1 局部場電位記錄與預處理
共有9名帕金森病患者參與本實驗，在英國牛津大學John Radcliffe Hospital完成丘腦底核深部腦刺激電極植入手術，植入電極為美敦力DBS 3389，電極直徑為1.27 mm，各電極觸點長度為1.5 mm，4個觸點之間的相鄰間隔為0.5 mm，電極植入后通過核磁共振成像對位置進行確認。
電極植入4~6 d后，從外置的電極導線采用雙極方法(相鄰觸點0-1,1-2,2-3)記錄丘腦底核深部局部場電位。從震顫相關的前臂伸肌和屈肌肌群同步記錄表面肌電圖。信號采樣頻率均為2 000 Hz，信號經由CED1902隔離放大器放大(Cambridge Eelectronic Design，UK)，局部場電位放大10 000倍，肌電信號放大1 000倍，帶通濾波參數為0.5~500 Hz，并使用CED1401 II和Spike2軟件(Cambridge Eelectronic Design，UK)對信號進行模數轉換和采集。實驗獲得當地道德倫理委員會批準許可。數據已得到英國牛津大學功能神經外科研究組授權使用。
預處理對信號進行90 Hz抗混疊低通濾波(切比雪夫I型濾波器)，降采樣2 000 Hz信號至250 Hz，采用自適應迭代濾波除去50 Hz工頻干擾，最后進行3 Hz高通濾波去除低頻成分(切比雪夫I型濾波器)。
1.2 時頻分析
應用STFT方法觀察局部場電位不同成分的能量密度隨時間的變化情況。離散隨機信號x(i),i=1,2,3…的STFT定義為
| ${\rm{STFT}}\left( {k,m} \right) = \sum\limits_{i = 1}^L {x\left( i \right){h^ * }\left( {i - k} \right)W_L^{mi}} ,$ |
其中窗函數h(k)時刻k為中心，每一個時刻k都可以得到一個不同的頻譜，則時頻譜密度P(k,ω)為
| $P\left( {k,\omega } \right) = {\left| {{\rm{STFT}}\left( {k,m} \right)} \right|^2}$ |
1.3 跨頻率相關分析
根據不同患者信號的譜峰分布，采用帶通濾波提取震顫相關的頻率成分，通過Hilbert變換提取各成分的幅度包絡，進而用互相關函數定量分析各成分包絡信號之間的相關程度。對于有限長的樣本，互相關函數的估計為
| ${\hat r_{xy}}\left( m \right) = \frac{1}{N}\sum\limits_{n = 0}^{N - 1 - \left| m \right|} {x\left( n \right)y\left( {n + m} \right)} \;\;0 \le \left| m \right| \le N - 1$ |
其中x(i)和y(i)的長度均為N，互相關函數序列長度為2N-1。
1.4 時域Granger因果分析
在Hilbert變換提取不同頻率成分的幅度包絡信息基礎上，Granger因果分析可對兩組時間序列之間的定向線性因果關系進行度量，衡量一個時間序列x(t)是否對另一個時間序列y(t)的發生有貢獻，進而估計兩個信號之間存在的因果影響的程度。
Granger因果分析簡單定義如下，若兩個廣義平穩時間序列x(t)和y(t)之間的線性關系可以用雙變量AR模型來表示，即
| $\begin{array}{l} x\left( t \right) = \sum\limits_{k = 1}^P {{a_{11}}\left( k \right)x\left( {t - k} \right)} + \sum\limits_{k = 1}^P {{a_{12}}\left( k \right)y\left( {t - k} \right)} + {v_1}\left( t \right)\\ y\left( t \right) = \sum\limits_{k = 1}^P {{a_{21}}\left( k \right)x\left( {t - k} \right)} + \sum\limits_{k = 1}^P {{a_{22}}\left( k \right)y\left( {t - k} \right)} + {v_2}\left( t \right), \end{array}$ |
其中p為模型階數，aij(k)為模型參數，v1(t)和v2(t)分別為x(t)和y(t)的預測誤差。當x(t)和y(t)分別僅用自身作預測，v1(t)和v2(t)的方差分別為Σ1和Σ2。當x(t)和y(t)利用自身和對方的前p個值作預測，此時v1(t)和v2(t)的方差分別為Σ11和Σ22。對于y(t)引起x(t)的因果關系時域上可以定義為
| ${F_{y \to x}} = \ln \frac{{\sum 1 }}{{\sum 1 1}}\left( {{f_{y \to x}} \ge 0} \right)$ |
若Σ11<Σ1，則說明加入y(t)對預測x(t)有幫助，預測誤差減少，認為y(t)對x(t)的發生有貢獻，反之亦然。
1.5 雙譜分析
高階譜分析中的雙譜分析定義為隨機時間信號三階累積量的二維傅里葉變換：
| $\begin{array}{l} R\left( {{m_1},{m_2}} \right) = E\left[ {x\left( n \right)x\left( {n + {m_1}} \right)x\left( {n + {m_2}} \right)} \right]\\ B\left( {{\omega _1},{\omega _2}} \right) = \sum\limits_{{m_1} = - \infty }^\infty {\sum\limits_{{m_2} = - \infty }^\infty {R\left( {{m_1},{m_2}} \right)} } \exp \left[ { - j\left( {{\omega _1}{m_1} + {\omega _2}{m_2}} \right)} \right], \end{array}$ |
其中R(m1,m2)為隨機信號x(t)的三階累積量，m1和m2為時間延遲量，B(ω1,ω2)為頻率對(ω1,ω2)的雙譜值。
雙譜保留了信號的相位信息，可以檢測出兩個頻率成分與它們的和頻率成分之間有無二次相位耦合作用[13]。二次相位耦合指一個非線性作用過程中兩個頻率成分在同時滿足頻率耦合和相位耦合的條件下，會在它們的和頻率處激發產生能量的現象。由下式中兩個互相獨立的隨機信號y(t)和z(t)為例說明：
| $\begin{array}{l} y\left( t \right) = \cos \left( {2\pi {\omega _1}t + {\varphi _1}} \right) + \cos \left( {2\pi {\omega _2} + {\varphi _2}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\cos \left( {2\pi {\omega _3} + {\varphi _3}} \right)\\ z\left( t \right) = \cos \left( {2\pi {\omega _1}t + {\varphi _1}} \right) + \cos \left( {2\pi {\omega _2} + {\varphi _2}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\cos \left( {2\pi {\omega _3} + {\varphi _4}} \right), \end{array}$ |
y(t)和z(t)均滿足ω1+ω2=ω3,φ1,φ2,φ3是[-π,π]內均勻分布的獨立隨機變量，且φ1+φ2=φ4。信號z(t)中存在二次相位耦合作用，ω1,ω2兩個成分的二次相位耦合對ω3處能量的產生有貢獻作用，而信號y(t)中不存在二次相位耦合作用，ω3處能量產生與ω1,ω2成分無關。
采用間接法估計雙譜，選用Parzen窗，每段數據長2 s，時間延遲為1.2 s，數據重疊1 s，1 024點FFT。在實際應用中，受到估計誤差和信號長度的限制，雙譜域中可能出現僅有頻率耦合沒有相位耦合的偽峰。替代數據方法消除了原信號中原有的非線性過程，僅保留了與原信號相同的平均功率譜密度特征，是鑒別偽峰的有效方法[14]。采用AAFT方法產生替代數據[15]，將原數據產生的雙譜與替代數據產生的雙譜進行比較，若兩者之間有顯著性差異(P<0.05，t-test)，則代表有統計學意義的雙譜譜峰存在。
2 結果
2.1 譜分析與時頻分析
在深部腦刺激電極植入后，在震顫和靜息狀態下同步記錄肌電信號和丘腦底核局部場電位，功率譜分析顯示9名患者丘腦底核場電位信號在震顫狀態下均存在一定程度的似諧波現象，如表 1所示，所有患者均存在震顫頻率(基頻)的譜峰(3~5 Hz)和2倍頻譜峰，5名患者中存在3倍頻譜峰。
表1
9例帕金森病患者震顫狀態丘腦底核局部場電位譜峰分布
Table1.
Distribution of power spectral peaks in tremor LFPs recorded from STNs of 9 Parkinson’s disease patients
以#1的40 s的肌電信號和局部場電位信號，如圖 1(a)、(b)所示為例說明，0~28 s時間段帕金森病患者處于震顫狀態，28~40 s時間段帕金森病患者處于靜息狀態。震顫狀態下3.5、7.0、10.5 Hz處有明顯的功率譜峰，如圖 1(c)所示；靜息狀態沒有這種現象，如見圖 1(d)所示。圖 1(e)為#1局部場電位信號的STFT結果，在震顫狀態時，4個頻率成分均保持較高的能量密度，28 s后從震顫狀態轉為靜息狀態，10.5 Hz和14.0 Hz能量密度急劇下降，而3.5 Hz和7.0 Hz處能量密度則在32 s后出現減弱，在其余患者均可見震顫狀態下各倍頻成分能量密度高于靜息狀態，以及各譜峰幅度隨時間而動態變化的現象，進而應用相關分析和因果分析度量各譜峰幅度動態變化之間的關系。
圖1
#1局部場電位時頻分析
(a)震顫狀態(0~28 s)和靜息狀態(28~40 s)的相關肌電信號；(b)與A同步記錄的震顫-靜息狀態的丘腦底核局部場電位；(c)震顫狀態丘腦底核局部場電位平均功率譜密度；(d)靜息狀態丘腦底核局部場電位平均功率譜密度；(e) 是(b)中局部場電位的時頻譜圖
Figure1. Time-frequency analysis of LFP signals in #1 patientTremor (0~28 s) and resting (28~40 s) related EMGs (a),LFPs (b),power spectrum of tremor (c) and resting (d) states,short-time Fourier transform of LFPs (e)
2.2 局部場電位各成分相關分析
對9例帕金森病患者震顫狀態的場電位信號進行跨頻率相關分析。根據每例患者局部場電位信號譜峰分布位置不同，分別對各例患者信號進行帶通濾波，提取各基頻、2倍頻和3倍頻頻段成分。應用Hilbert變換提取各成分的幅度包絡，計算兩兩間的互相關函數。圖 2為9名患者各頻率成分包絡間的最大互相關系數，倍頻成分的包絡信號互相關系數在0.18~0.50范圍內，說明各倍頻成分的幅度包絡信號間保持一定的相關性。
圖2
震顫狀態下9例帕金森病患者各成分包絡之間的平均最大互相關系數(均值±標準差)
Figure2.
Average of maximum cross-correlation coefficients between different components in LFPs of 9 Parkinson’s disease patients during tremor state (mean±standard deviation)
2.3 局部場電位各成分Granger因果分析
應用帶通濾波和Hilbert變換得到場電位各頻率成分的包絡，進而可利用時域Granger因果分析的方法估計各成分間的因果關系，平均因果分析結果如圖 3所示。9例帕金森病患者信號的實驗結果平均值顯示，震顫基頻對2倍頻(顯著性檢驗P<0.05，paired t-test)、震顫基頻對3倍頻、2倍頻對3倍頻有較強的因果值，說明處于較高頻段的倍頻成分的產生受到處于低頻段的基頻或倍頻成分的影響較大。
圖3
9名帕金森病患者震顫狀態下場電位各成分包絡之間的Granger因果分析(均值±標準差)
Figure3.
Granger causality between different components in LFPs of 9 Parkinson’s disease patients during tremor state (mean±standard deviation)
2.4 局部場電位雙譜分析
對震顫和靜息狀態下的局部場電位進行雙譜分析和統計學檢驗。9例帕金森病患者局部場電位信號在震顫狀態下的二次相位耦合作用主要集中基頻與基頻、2倍頻、3倍頻之間，如圖 4所示，說明在震顫狀態下震顫頻率自身耦合對2倍頻處的能量峰值產生有貢獻作用，震顫頻率和2倍頻之間的二次相位耦合對3倍頻處的能量產生有貢獻作用。這證明了震顫頻率相關的神經波動的存在，以及與其他倍頻成分的耦合作用驅動了更高頻率，如beta頻段成分的產生。靜息狀態下非線性耦合作用在雙譜域上的分布位置則不能解釋倍頻成分之間的驅動關系。
圖4
9例患者震顫狀態(實心圓)和靜息狀態(空心圓)信號雙譜譜峰位置分布
Figure4.
Bispectral peak positions in tremor state (filled circles) and resting state (unfilled circles) of 9 Parkinson’s disease patients
以患者#1為例進一步說明。與替代數據結果比較進行統計學顯著性檢測，患者#1震顫狀態下的場電位雙譜分析結果表明有6個譜峰如圖 5(a)、(c)所示，分別位于雙譜域中基頻-基頻、基頻-2倍頻、基頻-3倍頻的相應位置(P<0.05)，而在靜息狀態下僅在基頻-基頻處有一相對減弱的譜峰(P<0.05)。這說明在震顫狀態下，震顫基頻的自身耦合對2倍頻處的能量峰值產生有貢獻作用，震顫基頻和倍頻之間的二次相位耦合對3倍頻處的能量產生有貢獻作用。
圖5
局部場電位雙譜分析結果
#1震顫狀態丘腦底核局部場電位的雙譜幅度(a)和譜峰位置(c)；#1靜息狀態丘腦底核局部場電位的雙譜幅度(b)和譜峰位置(d)
Figure5. Bispectral analysis of LFP signalsAmplitude (a) and peak positions (c) of bispectra in tremor state of #1; Amplitude (b) and peak positions (d) of bispectra in resting state of #1
3 討論與結論
本文研究了帕金森病患者丘腦底核局部場電位在震顫和靜息兩種帕金森病疾病狀態下神經波動成分之間的作用關系，在患者震顫狀態的局部場電位功率譜中發現了“m×震顫頻率”的似諧波現象(m=1,2,3)。進而應用跨頻率相關分析、時域Granger因果分析和雙譜分析，研究了各倍頻成分之間的線性和非線性作用關系。
時頻分析結果顯示，震顫基頻和2倍頻成分集中了震顫狀態主要的能量分布，基頻、2倍頻成分和3倍頻保持一定的同步性。跨頻率相關分析發現，各倍頻成分的幅度間保持0.18~0.50的相關性，表明各成分間存在相關關系，一定程度上解釋了信號時頻譜圖中能量密度變化的現象。時域Granger因果分析量化了一個成分對產生另一個成分的影響，
實驗結果顯示了低頻處的基頻和2倍頻成分對高頻處的3倍頻成分相比高頻成分對低頻成分有更強的線性因果關系。
雙譜分析檢測了倍頻成分間的二次非線性關系，解釋了震顫基頻自身、震顫基頻與2倍頻頻段倍頻成分因為相位耦合作用激發了2倍頻、3倍頻成分的能量，進而解釋了時頻分析中震顫到靜息狀態轉變過程中各成分能量密度的變化現象。
震顫是帕金森病的典型癥狀，其表現形式為有規律的肢體顫動，藥物治療對震顫的控制作用較為有限。研究人員對震顫的發病原因提出過多種假設，但至今其神經機制仍不清楚。本工作在我們以往的工作基礎上進一步證實震顫的產生與基底神經節神經活動的多個波動成分有關，而且這些成分并非獨立，存在著相互之間的線性和非線性交互作用，并具有隨時間改變的動態性，且與患者個體有關，而這些神經活動是由一個神經環路亦或是來自于多個神經環路匯集于基底神經節的仍待進一步探索。
