本文基于多重分形去趨勢波動分析方法（MFDFA）實現了對受試者心理壓力狀態的評估。研究針對不同心理壓力狀態下，腦電信號多重分形去趨勢波動分析中最優分形階數確定問題，重點分析了多重分形去趨勢波動分析方法中的奇異指數、Hurst 指數等參數與階數的關系，進而確定最優分形階數，實現了基于腦電信號多重分形去趨勢波動分析的心理壓力狀態評估。試驗采集了 14 名在校學生有/無心理壓力狀態下的腦電信號，分別比較了奇異指數、奇異維數、Hurst 指數、質量指數與階數關系，確定了最優分形階數范圍為 [—5, 5]，實現了基于腦電信號 β 波多重分形去趨勢波動分析方法的心理壓力狀態評估。研究結果表明，心理壓力狀態下，腦電信號的 Hurst 指數和質量指數大于無壓力狀態下，腦電信號的相應參數，隨著階數的變大，Hurst 指數減小，趨近于定值，而質量指數增大，奇異值隨階數的變化幅度較明顯。本文還比較了有/無心理壓力狀態下，腦電信號的峰值和奇異譜寬度，結果表明，不同心理壓力狀態下腦電信號多重分形譜特性不同，心理壓力狀態下，腦電信號的奇異譜寬度明顯大于無壓力狀態下腦電信號的奇異譜寬度。本文研究結果說明，該方法可以有效地評估心理壓力狀態，為實現心理壓力狀態干預，提高心理健康等提供支持與幫助。
引用本文： 李昕, 孫小棋, 齊曉英, 侯永捷, 田彥秀. 面向心理壓力評估的腦電信號多重分形去趨勢波動分析方法研究. 生物醫學工程學雜志, 2017, 34(2): 180-187. doi: 10.7507/1001-5515.201512028 復制
引言
腦電（electroencephalogram, EEG）信號是生物電活動和腦功能狀態的綜合反映，具有微弱、非線性和強噪聲的特點。由于腦電信號中蘊含有豐富的生理、心理和病理信息，如何充分挖掘這些信息與有效地提取腦電信號特征，成為腦電信號研究的關鍵問題。
目前常用的腦電信號特征提取方法包括：① 時頻分析。例如：小波變換和多分辨率特征分析方法。② 能量分析。例如：復雜度分析可以測量時間序列的規律性，反映腦電信號時間序列復雜程度；近似熵可以衡量時間序列的復雜性，時間序列越復雜則近似熵值就越大。而小波熵、功率譜熵等則主要研究腦電信號多頻率成分的混亂程度以及信號能量在頻率成分上的分布。③ 多重分形。多重分形去趨勢波動分析方法（multi-fractal de-trended fluctuation analysis, MFDFA）主要分析腦電信號分形結構，挖掘腦電信號中長程相關性和冪指數規律等特性。其中，多重分形方法是用一個譜函數從不同層次說明分形體的生長特征，并反映大腦接受或處理信息的有序程度，是分析腦電信號特征的一種有效方法。
1967 年，Mandelbrot 在美國《科學》雜志上發表關于“英國的海岸線有多長”的一篇論文，提出分形概念。1986 年 Halsey 等[1]提出多重分形概念，并指出腦電信號具有長程相關性和冪律指數規律，即腦電信號具有標度不變性。1994 年 Peng 等[2]提出去趨勢波動分析（de-trended fluctuation analysis, DFA）方法。1998 年，Watters 采用 DFA 方法檢測腦電信號的相關性，結果顯示人閉目靜坐時腦電信號的標度指數為 1.26，表明腦電信號具有分形特性。2002 年 Kantelhardt 等[3]在 DFA 的基礎上提出多重分形去趨勢波動分析方法，發現非平穩有限時間序列不僅顯示了分形結構特征的長程相關性和標度不變性，還具有多重分形屬性。這種方法是 DFA 方法的一種直接推廣，是優于標準多重分形形式體系的多重分析方法，可以運用于帶有趨勢的非平穩時間序列的分析。
鄒鳴等[4]利用多重分形去趨勢波動分析方法研究睜眼、閉眼、正常人、癲癇未發作期和發作期患者的腦電信號變化，通過比較奇異譜寬度和 Hurst 指數區分不同的腦電信號狀態，得到了較好的分析結果。Bachmann 等[5]利用 DFA 方法對沮喪狀態下采集的腦電信號進行分析，結果表明沮喪狀態下的腦電信號具有多重分形特性，而且相對于正常人，其腦電信號時間序列具有持久的長程相關性，心情沮喪時的腦電信號時間序列具有相對短暫的長程相關性。Dutta 等[6]討論了多重分形參數是否可以作為健康人和癲癇患者的區分指標，研究表明，癲癇患者發病時，腦電信號的多重分形程度比正常人高出很多，并且正常人在睜眼和閉眼時的多重分形程度也不相同，證實了利用多重分形參數在癲癇診斷中應用的可行性。Marton 等[7]將 DFA 方法應用到腦電信號的研究，比較了分形結構在腦電信號低頻成分和高頻成分中的不同。
多重分形去趨勢波動分析可以有效提取腦電信號分形特征。目前，該方法應用在心理壓力評估中的研究還不多見。本文基于多重分形去趨勢波動分析方法，分析心理壓力狀態下的腦電信號，在分析確定最優分形階數的基礎上，提取腦電信號分形特征，實現心理壓力狀態的評估。該方法可以有效地評估心理壓力狀態，為進一步實現心理壓力狀態干預，提高心理健康提供支持與幫助。
1 多重分形去趨勢波動分析方法
多重分形去趨勢波動分析方法是 Kante-lhardt[3]等在 DFA 方法的基礎上提出的一種新的多重分形分析方法。對長度為n 非平穩時間序列 {xk }，k=1, 2, 3, ,n 計算分形指數，步驟如下：
（1）計算該時間序列的均值：
| $\bar x = \frac{1}{n}\sum\limits_k^n {{x_k}} $ |
（2）計算累計離差，即時間序列減去均值：
| $Y(i) = \sum\limits_{k = 1}^i {[{x_k} - \overline x ]} \;\;\;\;\;i = 1,{\rm{ }}2,{\rm{ }}3,{\rm{ }}\cdot \! \cdot \! \cdot{\rm{ }},n$ |
基于以上處理，噪聲時間序列轉化為隨機游走時間序列 {Y（i）}。進一步對該時間序列進行多重分形分析。
（3）將時間序列 {Y（i）}分成m 個長度為s 的不重疊區間，每一個區間v 含有s 個點，對每一個區間采用最小二乘法進行多項式擬合，得到n階擬合多項式，稱為擬合趨勢值：
| $\begin{aligned}& {y_v}(i) = {a_1}{i^k} + {a_2}{i^k} + \cdot \! \cdot \! \cdot{a_k}i + {a_{k + 1}} \\ &i = 1,{\rm{ }}2,{\rm{ }}3,{\rm{ }} \cdot \! \cdot \cdot ,s,k = 1,{\rm{ }}2,{\rm{ }} \cdot \! \cdot \! \cdot {\rm{ }},n\end{aligned}$ |
（4）每一個區間v 上所有點的值減去擬合趨勢值yv （i），得到擬合殘差，即區間的局部趨勢，該趨勢反映變化幅度敏感，而不是變化結構。計算局部趨勢的均方值得到尺度函數F（v,s）。
| ${F^2}(v,s) = \frac{1}{s}\sum\limits_{i = 1}^s {\{ Y[(v - 1)s + i]} - {y_v}(i){\} ^2}$ |
（5）由于區間中含有的點數不同，若時間段較短，即所含數據點較少，則快波影響尺度函數，反之，則慢波影響尺度函數，所以在計算尺度函數時，將時間序列分為多段不同長度的尺寸，突出時間序列中快波和慢波的影響，記為RMS{ns}（v）。
| $RMS\{ ns\} (v) = F(v,s)$ |
其中ns 代表時間段不同的個數，v 代表一個區間。
（6）計算所有區間的平均值得到q 階波動函數。
| ${F_q}(s) = {\{ \frac{1}{{{n_s}}}\sum\limits_{v = 1}^{{N_s}} {[{F^2}} (v,s)]^{q/2}}{\} ^{1/q}}$ |
一般情況q 變量的取值可以是任意值，例如q=2 時相當于標準的 DFA 程序，波動函數隨著s 的增加而增加，一般s≥m。
（7）若時間序列 {xk } 是長程相關的，則：
| ${F_q}(s) \sim {s^{h(q)}}$ |
對于一些大的尺度范圍，例如s>n/4，Fq （s）會變得不穩定。
1.1 Hurst 指數
Hurst 指數是英國 H.E Hurst 提出的一種判別時間序列是否對時間具有依賴性的參數，主要用來判斷數據是否具有混沌性以及表現序列特征的微小變化[9]。根據 Hurst 指數值范圍不同，數據相關性也不同，主要分為 3 個不同類型：
（1）0<Hurst<0.5，表示時間序列具有長期相關性，但總體趨勢與以前相反，即反持續，這種反持續依賴于 Hurst 值與 0 的距離，越接近 0 表示越負相關。
（2）Hurst=0.5，表示時間序列是不相關的。
（3）Hurst>0.5，時間序列長期相關，并且具有持續性，Hurst 值越接近 1，相關性越強。
在多重分形去趨勢波動分析方法中，由式（7）可以看出波動函數與每個區間所含點數是指數關系，兩邊取對數：
| $\lg {F_q}(s) = h(q)\lg s + C$ |
| $h(q) \sim q$ |
由式（8）計算線性回歸得到的斜率為標度指數h（q），即 Hurst 指數。根據q 取值可以描述不同幅度波動的標度行為，例如q 取較大值時，描述幅度波動較大時的標度行為，相對于較小的q 值描述幅度波動較小的標度行為。單一階數對應一個 Hurst 指數，依賴于q，稱為q-Hurst 指數，如式（9）所示。例如q=2時，標度指數h（2）稱為經典 Hurst 指數。
1.2 質量指數
質量指數τ（q），即多重分形過程中的尺度指數之一，是分形行為的另一個重要刻畫指標。多重分形去趨勢波動分析方法第一個經典過程是將q 階h（q）轉化為q 階質量指數τ（q），如式（10）所示：
| ${τ} (q) = qh(q) - 1$ |
如式（10）所示，質量指數與階數和 Hurst 指數相關，同時說明了q-Hurst 指數和質量指數這兩種多重分形尺度指數之間的關系。質量指數一般用于計算奇異指數和奇異維數。質量指數是非線性函數，若質量指數為線性函數，則不是多重分形，而是單分形過程。
1.3 分形奇異維數與奇異指數
奇異指數α 用來描述區間不同的奇異程度，與奇異性成反比。分形奇異維數f（α）是建立在奇異強度的基礎上，直接反應了奇異指數的分形維數，表示相應奇異指數的分布密度。
奇異指數α 和q 階奇異維數f（α），如式（11）所示：
| $\begin{array}{l}{{α}} = {τ} '(q)\\f({{α}}) = q{{α}} - {τ} (q)\end{array}$ |
α 反映各個區間不同的奇異程度，f（α）反映了具有奇異指數α 的分形維數。α–f（α）曲線一般稱為多重分形譜曲線，形狀為單峰拱形。
△α=αmax–αmin 表示分形譜的寬度，代表最大概率和最小概率之間的差別，△α 越大，時間序列分布越不均勻，多重分形強度越大。
2 數據處理
2.1 數據采集
受試者為 14 名在校大學生（7 男 7 女），身體健康，且均為右利手。試驗在一個隔絕聲音、隔絕光線的屏蔽室內進行，試驗之前受試者頭發洗凈、吹干并保持良好狀態；試驗控制在進餐 3 小時內進行，以減少低血糖的影響。限時任務為背誦古文材料。受試者靜息 5 min，腦電信號波形平穩后，睜眼靜息 2 min 后開始進行限時背誦，時長為 1 min。共采集 92 組腦電信號數據，包括 44 組靜息樣本和 48 個游戲樣本。
腦電信號采集儀器使用的是美國生產的 NeuroScan 64 導腦電信號采集系統，電極按照 10～20 國際電極系統放置。信號采樣率設為 1 000 Hz，0～30 Hz 低通濾波。參考電極是左右耳乳突，使用儀器本身自帶的分析軟件去除眼動偽跡，使用 eeglab 工具箱進行 128 Hz 重采樣，使用獨立分量分析（independent component analysis, ICA）方法對數據進行去噪處理。截取不同狀態下的平穩腦電信號數據進行下一步分析。
2.2 腦電信號多重分形去趨勢波動參數分析
2.2.1 多重分形去趨勢波動參數 試驗設計是速記試驗，受試者為成年人。試驗過程中受試者處于睜眼狀態，α 波會減弱或者消失，不宜用于試驗分析。對于 θ 波和 δ 波來說，嬰兒時期這兩種波占據主要地位，而成年人這兩種節律波只有在困倦或者睡眠狀態下才會出現。本文試驗是在睜眼狀態下進行的，所以 θ 波和 δ 波不分析。大腦皮層處于興奮狀態時 β 波變化明顯，所以主要提取腦電信號中的 β 波用于壓力評估與分析。
首先，分析有/無壓力腦電信號的分形特征，如圖 1 所示的多重分形譜圖可以看出，有/無壓力狀態下腦電信號多重分形譜曲線形狀均為單峰拱形，說明有/無壓力腦電信號具有多重分形特性。由多重分形譜曲線可知有壓力腦電信號的奇異譜寬度大于無壓力腦電信號的奇異譜寬度。
圖1
多重分形階數與多重分形指數的關系
Figure1.
Relationship between multifractal spectrum orders and multifractal spectrum indexes
進一步，為了確定最優分形階數，分析了 Hurst 指數、質量指數與分形階數的關系以及奇異指數、奇異維數與階數關系。如圖 1 所示，實線代表限時背誦腦電信號數據，即有壓力腦電信號數據，虛線表示靜息腦電信號數據。有壓力腦電信號的 Hurst 指數大于無壓力腦電信號。隨著階數的增大，Hurst 指數兩條線趨于定值，且大于 0.5，并且階數接近于 0 時，Hurst 接近于 1，腦電信號相關性最強；有壓力腦電信號質量指數曲率大于無壓力腦電信號，說明有壓力腦電信號分形特性較明顯。
奇異指數和分形維數隨階數的變化規律亦如圖 1 所示，兩者與階數都是非線性關系，證明利用多重分形去趨勢波動方法分析壓力腦電信號是可行的。
2.2.2 基于多重分形去趨勢波動參數的階數選擇
階數選擇是多重分形的關鍵，對于試驗結果有著較大的影響。階數選擇過大或過小都有可能造成結果的不準確性，影響數據分析的效果。作為判斷數據混沌性的一個指標，首先分析 Hurst 指數與不同階數之間的變化關系。
Hurst 指數與階數關系如圖 2 所示，實線代表限時背誦，虛線代表靜息狀態。Hurst 指數隨著階數增大呈減小趨勢，凸顯出在階數大于 10 或小于 —10 時有/無壓力腦電信號的 Hurst 指數幾乎不變，接近常數，并且在大于 10 范圍內兩者 Hurst 指數接近相等。而在 [0, 8] 范圍內，Hurst 指數接近 1，且在大于 0.5 范圍內，所以針對 Hurst 指數階數選擇在 [—5, 8] 之間。
圖2
Hurst 指數與階數關系
Figure2.
Relationship between different orders and Hurst index
質量指數與階數關系如圖 3 所示，曲線變化趨勢基本不變，但質量指數初值隨階數增大在逐漸減小，如式（11）所示，初值影響了奇異維數的變化。有壓力腦電信號曲線曲率較無壓力腦電信號略大，并且兩種狀態的質量指數隨著階數增大趨于相等，所以針對質量指數，階數不宜過大。由式（11）可知，質量指數曲線的曲率與奇異指數和奇異維數有關，所以階數的選擇依據階數與奇異指數和奇異維數的關系。
圖3
質量指數與階數關系
Figure3.
Relationship between different orders and quality index
奇異指數是用來描述復雜體系中各個區間不同的奇異程度，奇異指數越小，奇異性越大，奇異指數與階數關系如圖 4 所示。奇異指數的初值隨階數增大幾乎不變，但在大于 5 和小于 —5 范圍內奇異指數隨階數的變化不明顯。例如：階數大于 5 時，有/無壓力的奇異指數較為接近，不容易區分，所以階數最好選在 [—5, 5]。
圖4
奇異指數與階數關系
Figure4.
Relationship between different orders and singular index
多重分形譜反映奇異指數與奇異維數之間關系，表示奇異值所在集合的差別，奇異維數與階數的關系如圖 5 所示，隨著階數增大，奇異維數趨于定值。通過觀察圖 5 發現有/無壓力的奇異維數在適當的階數范圍內區分較明顯，所以取中間階數 [—8, 8] 適宜。
圖5
奇異維數與階數關系
Figure5.
Relationship between different orders and singular dimension
不同階數奇異維數與奇異指數關系變化如圖 6 所示，奇異指數和奇異維數關系幾乎不隨階數范圍的增大而有所變化。比較圖 6 中三個圖發現階數在 [—10, 10]，奇異譜寬度更容易看出。
圖6
奇異維數與奇異指數關系圖
Figure6.
Relationship between singular index and singular dimension
綜上，質量指數和奇異維數隨階數變化影響不明顯，主要為 Hurst 指數和奇異指數，綜合以上分析結果，選取階數范圍為 [—5, 5]。
3 基于多重分形去趨勢波動分析方法心理壓力評估
通過心理壓力腦電信號多重分形去趨勢波動分析參數分析，確定最優階數范圍為[—5, 5]。在這個范圍內不僅 Hurst 指數趨近于 1，說明有/無壓力狀態腦電信號均為長相關關系，而且奇異譜曲線都為單峰曲線，說明有/無壓力腦電信號具有多重分形特性，基于以上分析，本文采用奇異譜寬度評估心理壓力狀態。
3.1 算法流程
算法總流程如圖 7 所示
圖7
算法流程圖
Figure7.
Flow chart of multifractal de-trended fluctuation analysis
3.2 結果與分析
首先對采集到的腦電信號數據進行篩選。分形是對具有自相似性無固定特征長度圖形結構的總稱，簡單來說分形可以理解為一座橋：一岸是單個神經元行為，另一岸代表整個大腦的行為。腦電信號的自相似性來源于其本身分形特性，由文獻[8]可知，正常人清醒腦電信號多重分形譜圖為單峰拱形。
如式（10）、式（11）所示，可以看出一個階數對應唯一一個奇異指數和奇異維數，具有唯一性。正常人的奇異譜曲線為單峰拱形，通過觀察采集數據的奇異譜曲線效果對數據進行了篩選，不正常奇異譜曲線如圖 8 所示。數據多重分形奇異譜曲線若出現類似圖 8 的曲線，則刪去該組數據。最終選取 22 組有/無壓力腦電信號數據進行分析。同時對選取的數據進行了統計分析，如表 1、表 2 所示。腦電信號多重分形譜寬度分析結果如表 3 所示。
圖8
不正常奇異譜曲線
Figure8.
Abnormal singular spectrum curve diagram
如表 1 所示，是兩組數據的統計特征不同。如表 2 所示，分析了 Levene 方差齊性檢驗，檢驗兩組方差是否齊，F=0.55，P=0.46 說明方差是齊的。均值方程t 檢驗中P=0.004，說明結果具有統計學意義。
表1
腦電信號統計分析
Table1.
Statistics of two groups of EEG
表2
獨立樣本t 檢驗
Table2.
Independent samplet test
表3
腦電信號多重分形奇異指數
Table3.
Singular index of multi-fractal de-trended fuctuation analysis with electroencephalogram
不同的人對壓力反應不同，所表現出來的多重分形譜不同，即存在個體差異。分析過程中，減去偏差較大的數據，選取 22 組數據進行分析。分別求有/無壓力腦電信號數據奇異譜寬度的最大值、最小值以及平均值并進行比較，比較結果如表 4 所示。
表4
奇異譜寬度比較
Table4.
Comparison of singular spectrum width
有壓力腦電信號的奇異譜寬度最大值、最小值和平均值分別為 0.81、0.35 和 0.59，均高于無壓力腦電信號的奇異譜寬度。基于以上分析，在有效分析得到最優分形階數的基礎上，基于多重分形譜寬度等參數能夠有效分析有/無壓力腦電信號數據，進而有效評估有/無心理壓力狀態。
4 結論
本文基于多重分形去趨勢分析方法實現腦電信號分析與壓力狀態評估。分析了多重分形去趨勢波動分析方法中的質量指數、奇異維數、奇異指數、Hurst 指數與階數的關系。分析結果表明：質量指數和奇異維數隨階數變化不明顯，而 Hurst 指數和奇異指數隨階數變化明顯，心理壓力腦電信號分形階數選擇為 [—5, 5] 時，能夠充分挖掘腦電信號特征。
有壓力腦電信號的奇異譜寬度的最大值和最小值均明顯高于無壓力腦電信號的奇異譜寬度。奇異譜寬度平均值結果顯示，有壓力狀態下，奇異譜寬度平均值明顯高于無壓力狀態下奇異譜寬度平均值，基于奇異譜寬度分析，可以確定壓力狀態。基于多重分形去趨勢波動分析方法，可以有效地識別心理壓力狀態，進而通過相應的干預和預防能夠幫助改善心理壓力狀態，并創建健康的生活方式。
引言
腦電（electroencephalogram, EEG）信號是生物電活動和腦功能狀態的綜合反映，具有微弱、非線性和強噪聲的特點。由于腦電信號中蘊含有豐富的生理、心理和病理信息，如何充分挖掘這些信息與有效地提取腦電信號特征，成為腦電信號研究的關鍵問題。
目前常用的腦電信號特征提取方法包括：① 時頻分析。例如：小波變換和多分辨率特征分析方法。② 能量分析。例如：復雜度分析可以測量時間序列的規律性，反映腦電信號時間序列復雜程度；近似熵可以衡量時間序列的復雜性，時間序列越復雜則近似熵值就越大。而小波熵、功率譜熵等則主要研究腦電信號多頻率成分的混亂程度以及信號能量在頻率成分上的分布。③ 多重分形。多重分形去趨勢波動分析方法（multi-fractal de-trended fluctuation analysis, MFDFA）主要分析腦電信號分形結構，挖掘腦電信號中長程相關性和冪指數規律等特性。其中，多重分形方法是用一個譜函數從不同層次說明分形體的生長特征，并反映大腦接受或處理信息的有序程度，是分析腦電信號特征的一種有效方法。
1967 年，Mandelbrot 在美國《科學》雜志上發表關于“英國的海岸線有多長”的一篇論文，提出分形概念。1986 年 Halsey 等[1]提出多重分形概念，并指出腦電信號具有長程相關性和冪律指數規律，即腦電信號具有標度不變性。1994 年 Peng 等[2]提出去趨勢波動分析（de-trended fluctuation analysis, DFA）方法。1998 年，Watters 采用 DFA 方法檢測腦電信號的相關性，結果顯示人閉目靜坐時腦電信號的標度指數為 1.26，表明腦電信號具有分形特性。2002 年 Kantelhardt 等[3]在 DFA 的基礎上提出多重分形去趨勢波動分析方法，發現非平穩有限時間序列不僅顯示了分形結構特征的長程相關性和標度不變性，還具有多重分形屬性。這種方法是 DFA 方法的一種直接推廣，是優于標準多重分形形式體系的多重分析方法，可以運用于帶有趨勢的非平穩時間序列的分析。
鄒鳴等[4]利用多重分形去趨勢波動分析方法研究睜眼、閉眼、正常人、癲癇未發作期和發作期患者的腦電信號變化，通過比較奇異譜寬度和 Hurst 指數區分不同的腦電信號狀態，得到了較好的分析結果。Bachmann 等[5]利用 DFA 方法對沮喪狀態下采集的腦電信號進行分析，結果表明沮喪狀態下的腦電信號具有多重分形特性，而且相對于正常人，其腦電信號時間序列具有持久的長程相關性，心情沮喪時的腦電信號時間序列具有相對短暫的長程相關性。Dutta 等[6]討論了多重分形參數是否可以作為健康人和癲癇患者的區分指標，研究表明，癲癇患者發病時，腦電信號的多重分形程度比正常人高出很多，并且正常人在睜眼和閉眼時的多重分形程度也不相同，證實了利用多重分形參數在癲癇診斷中應用的可行性。Marton 等[7]將 DFA 方法應用到腦電信號的研究，比較了分形結構在腦電信號低頻成分和高頻成分中的不同。
多重分形去趨勢波動分析可以有效提取腦電信號分形特征。目前，該方法應用在心理壓力評估中的研究還不多見。本文基于多重分形去趨勢波動分析方法，分析心理壓力狀態下的腦電信號，在分析確定最優分形階數的基礎上，提取腦電信號分形特征，實現心理壓力狀態的評估。該方法可以有效地評估心理壓力狀態，為進一步實現心理壓力狀態干預，提高心理健康提供支持與幫助。
1 多重分形去趨勢波動分析方法
多重分形去趨勢波動分析方法是 Kante-lhardt[3]等在 DFA 方法的基礎上提出的一種新的多重分形分析方法。對長度為n 非平穩時間序列 {xk }，k=1, 2, 3, ,n 計算分形指數，步驟如下：
（1）計算該時間序列的均值：
| $\bar x = \frac{1}{n}\sum\limits_k^n {{x_k}} $ |
（2）計算累計離差，即時間序列減去均值：
| $Y(i) = \sum\limits_{k = 1}^i {[{x_k} - \overline x ]} \;\;\;\;\;i = 1,{\rm{ }}2,{\rm{ }}3,{\rm{ }}\cdot \! \cdot \! \cdot{\rm{ }},n$ |
基于以上處理，噪聲時間序列轉化為隨機游走時間序列 {Y（i）}。進一步對該時間序列進行多重分形分析。
（3）將時間序列 {Y（i）}分成m 個長度為s 的不重疊區間，每一個區間v 含有s 個點，對每一個區間采用最小二乘法進行多項式擬合，得到n階擬合多項式，稱為擬合趨勢值：
| $\begin{aligned}& {y_v}(i) = {a_1}{i^k} + {a_2}{i^k} + \cdot \! \cdot \! \cdot{a_k}i + {a_{k + 1}} \\ &i = 1,{\rm{ }}2,{\rm{ }}3,{\rm{ }} \cdot \! \cdot \cdot ,s,k = 1,{\rm{ }}2,{\rm{ }} \cdot \! \cdot \! \cdot {\rm{ }},n\end{aligned}$ |
（4）每一個區間v 上所有點的值減去擬合趨勢值yv （i），得到擬合殘差，即區間的局部趨勢，該趨勢反映變化幅度敏感，而不是變化結構。計算局部趨勢的均方值得到尺度函數F（v,s）。
| ${F^2}(v,s) = \frac{1}{s}\sum\limits_{i = 1}^s {\{ Y[(v - 1)s + i]} - {y_v}(i){\} ^2}$ |
（5）由于區間中含有的點數不同，若時間段較短，即所含數據點較少，則快波影響尺度函數，反之，則慢波影響尺度函數，所以在計算尺度函數時，將時間序列分為多段不同長度的尺寸，突出時間序列中快波和慢波的影響，記為RMS{ns}（v）。
| $RMS\{ ns\} (v) = F(v,s)$ |
其中ns 代表時間段不同的個數，v 代表一個區間。
（6）計算所有區間的平均值得到q 階波動函數。
| ${F_q}(s) = {\{ \frac{1}{{{n_s}}}\sum\limits_{v = 1}^{{N_s}} {[{F^2}} (v,s)]^{q/2}}{\} ^{1/q}}$ |
一般情況q 變量的取值可以是任意值，例如q=2 時相當于標準的 DFA 程序，波動函數隨著s 的增加而增加，一般s≥m。
（7）若時間序列 {xk } 是長程相關的，則：
| ${F_q}(s) \sim {s^{h(q)}}$ |
對于一些大的尺度范圍，例如s>n/4，Fq （s）會變得不穩定。
1.1 Hurst 指數
Hurst 指數是英國 H.E Hurst 提出的一種判別時間序列是否對時間具有依賴性的參數，主要用來判斷數據是否具有混沌性以及表現序列特征的微小變化[9]。根據 Hurst 指數值范圍不同，數據相關性也不同，主要分為 3 個不同類型：
（1）0<Hurst<0.5，表示時間序列具有長期相關性，但總體趨勢與以前相反，即反持續，這種反持續依賴于 Hurst 值與 0 的距離，越接近 0 表示越負相關。
（2）Hurst=0.5，表示時間序列是不相關的。
（3）Hurst>0.5，時間序列長期相關，并且具有持續性，Hurst 值越接近 1，相關性越強。
在多重分形去趨勢波動分析方法中，由式（7）可以看出波動函數與每個區間所含點數是指數關系，兩邊取對數：
| $\lg {F_q}(s) = h(q)\lg s + C$ |
| $h(q) \sim q$ |
由式（8）計算線性回歸得到的斜率為標度指數h（q），即 Hurst 指數。根據q 取值可以描述不同幅度波動的標度行為，例如q 取較大值時，描述幅度波動較大時的標度行為，相對于較小的q 值描述幅度波動較小的標度行為。單一階數對應一個 Hurst 指數，依賴于q，稱為q-Hurst 指數，如式（9）所示。例如q=2時，標度指數h（2）稱為經典 Hurst 指數。
1.2 質量指數
質量指數τ（q），即多重分形過程中的尺度指數之一，是分形行為的另一個重要刻畫指標。多重分形去趨勢波動分析方法第一個經典過程是將q 階h（q）轉化為q 階質量指數τ（q），如式（10）所示：
| ${τ} (q) = qh(q) - 1$ |
如式（10）所示，質量指數與階數和 Hurst 指數相關，同時說明了q-Hurst 指數和質量指數這兩種多重分形尺度指數之間的關系。質量指數一般用于計算奇異指數和奇異維數。質量指數是非線性函數，若質量指數為線性函數，則不是多重分形，而是單分形過程。
1.3 分形奇異維數與奇異指數
奇異指數α 用來描述區間不同的奇異程度，與奇異性成反比。分形奇異維數f（α）是建立在奇異強度的基礎上，直接反應了奇異指數的分形維數，表示相應奇異指數的分布密度。
奇異指數α 和q 階奇異維數f（α），如式（11）所示：
| $\begin{array}{l}{{α}} = {τ} '(q)\\f({{α}}) = q{{α}} - {τ} (q)\end{array}$ |
α 反映各個區間不同的奇異程度，f（α）反映了具有奇異指數α 的分形維數。α–f（α）曲線一般稱為多重分形譜曲線，形狀為單峰拱形。
△α=αmax–αmin 表示分形譜的寬度，代表最大概率和最小概率之間的差別，△α 越大，時間序列分布越不均勻，多重分形強度越大。
2 數據處理
2.1 數據采集
受試者為 14 名在校大學生（7 男 7 女），身體健康，且均為右利手。試驗在一個隔絕聲音、隔絕光線的屏蔽室內進行，試驗之前受試者頭發洗凈、吹干并保持良好狀態；試驗控制在進餐 3 小時內進行，以減少低血糖的影響。限時任務為背誦古文材料。受試者靜息 5 min，腦電信號波形平穩后，睜眼靜息 2 min 后開始進行限時背誦，時長為 1 min。共采集 92 組腦電信號數據，包括 44 組靜息樣本和 48 個游戲樣本。
腦電信號采集儀器使用的是美國生產的 NeuroScan 64 導腦電信號采集系統，電極按照 10～20 國際電極系統放置。信號采樣率設為 1 000 Hz，0～30 Hz 低通濾波。參考電極是左右耳乳突，使用儀器本身自帶的分析軟件去除眼動偽跡，使用 eeglab 工具箱進行 128 Hz 重采樣，使用獨立分量分析（independent component analysis, ICA）方法對數據進行去噪處理。截取不同狀態下的平穩腦電信號數據進行下一步分析。
2.2 腦電信號多重分形去趨勢波動參數分析
2.2.1 多重分形去趨勢波動參數 試驗設計是速記試驗，受試者為成年人。試驗過程中受試者處于睜眼狀態，α 波會減弱或者消失，不宜用于試驗分析。對于 θ 波和 δ 波來說，嬰兒時期這兩種波占據主要地位，而成年人這兩種節律波只有在困倦或者睡眠狀態下才會出現。本文試驗是在睜眼狀態下進行的，所以 θ 波和 δ 波不分析。大腦皮層處于興奮狀態時 β 波變化明顯，所以主要提取腦電信號中的 β 波用于壓力評估與分析。
首先，分析有/無壓力腦電信號的分形特征，如圖 1 所示的多重分形譜圖可以看出，有/無壓力狀態下腦電信號多重分形譜曲線形狀均為單峰拱形，說明有/無壓力腦電信號具有多重分形特性。由多重分形譜曲線可知有壓力腦電信號的奇異譜寬度大于無壓力腦電信號的奇異譜寬度。
圖1
多重分形階數與多重分形指數的關系
Figure1.
Relationship between multifractal spectrum orders and multifractal spectrum indexes
進一步，為了確定最優分形階數，分析了 Hurst 指數、質量指數與分形階數的關系以及奇異指數、奇異維數與階數關系。如圖 1 所示，實線代表限時背誦腦電信號數據，即有壓力腦電信號數據，虛線表示靜息腦電信號數據。有壓力腦電信號的 Hurst 指數大于無壓力腦電信號。隨著階數的增大，Hurst 指數兩條線趨于定值，且大于 0.5，并且階數接近于 0 時，Hurst 接近于 1，腦電信號相關性最強；有壓力腦電信號質量指數曲率大于無壓力腦電信號，說明有壓力腦電信號分形特性較明顯。
奇異指數和分形維數隨階數的變化規律亦如圖 1 所示，兩者與階數都是非線性關系，證明利用多重分形去趨勢波動方法分析壓力腦電信號是可行的。
2.2.2 基于多重分形去趨勢波動參數的階數選擇
階數選擇是多重分形的關鍵，對于試驗結果有著較大的影響。階數選擇過大或過小都有可能造成結果的不準確性，影響數據分析的效果。作為判斷數據混沌性的一個指標，首先分析 Hurst 指數與不同階數之間的變化關系。
Hurst 指數與階數關系如圖 2 所示，實線代表限時背誦，虛線代表靜息狀態。Hurst 指數隨著階數增大呈減小趨勢，凸顯出在階數大于 10 或小于 —10 時有/無壓力腦電信號的 Hurst 指數幾乎不變，接近常數，并且在大于 10 范圍內兩者 Hurst 指數接近相等。而在 [0, 8] 范圍內，Hurst 指數接近 1，且在大于 0.5 范圍內，所以針對 Hurst 指數階數選擇在 [—5, 8] 之間。
圖2
Hurst 指數與階數關系
Figure2.
Relationship between different orders and Hurst index
質量指數與階數關系如圖 3 所示，曲線變化趨勢基本不變，但質量指數初值隨階數增大在逐漸減小，如式（11）所示，初值影響了奇異維數的變化。有壓力腦電信號曲線曲率較無壓力腦電信號略大，并且兩種狀態的質量指數隨著階數增大趨于相等，所以針對質量指數，階數不宜過大。由式（11）可知，質量指數曲線的曲率與奇異指數和奇異維數有關，所以階數的選擇依據階數與奇異指數和奇異維數的關系。
圖3
質量指數與階數關系
Figure3.
Relationship between different orders and quality index
奇異指數是用來描述復雜體系中各個區間不同的奇異程度，奇異指數越小，奇異性越大，奇異指數與階數關系如圖 4 所示。奇異指數的初值隨階數增大幾乎不變，但在大于 5 和小于 —5 范圍內奇異指數隨階數的變化不明顯。例如：階數大于 5 時，有/無壓力的奇異指數較為接近，不容易區分，所以階數最好選在 [—5, 5]。
圖4
奇異指數與階數關系
Figure4.
Relationship between different orders and singular index
多重分形譜反映奇異指數與奇異維數之間關系，表示奇異值所在集合的差別，奇異維數與階數的關系如圖 5 所示，隨著階數增大，奇異維數趨于定值。通過觀察圖 5 發現有/無壓力的奇異維數在適當的階數范圍內區分較明顯，所以取中間階數 [—8, 8] 適宜。
圖5
奇異維數與階數關系
Figure5.
Relationship between different orders and singular dimension
不同階數奇異維數與奇異指數關系變化如圖 6 所示，奇異指數和奇異維數關系幾乎不隨階數范圍的增大而有所變化。比較圖 6 中三個圖發現階數在 [—10, 10]，奇異譜寬度更容易看出。
圖6
奇異維數與奇異指數關系圖
Figure6.
Relationship between singular index and singular dimension
綜上，質量指數和奇異維數隨階數變化影響不明顯，主要為 Hurst 指數和奇異指數，綜合以上分析結果，選取階數范圍為 [—5, 5]。
3 基于多重分形去趨勢波動分析方法心理壓力評估
通過心理壓力腦電信號多重分形去趨勢波動分析參數分析，確定最優階數范圍為[—5, 5]。在這個范圍內不僅 Hurst 指數趨近于 1，說明有/無壓力狀態腦電信號均為長相關關系，而且奇異譜曲線都為單峰曲線，說明有/無壓力腦電信號具有多重分形特性，基于以上分析，本文采用奇異譜寬度評估心理壓力狀態。
3.1 算法流程
算法總流程如圖 7 所示
圖7
算法流程圖
Figure7.
Flow chart of multifractal de-trended fluctuation analysis
3.2 結果與分析
首先對采集到的腦電信號數據進行篩選。分形是對具有自相似性無固定特征長度圖形結構的總稱，簡單來說分形可以理解為一座橋：一岸是單個神經元行為，另一岸代表整個大腦的行為。腦電信號的自相似性來源于其本身分形特性，由文獻[8]可知，正常人清醒腦電信號多重分形譜圖為單峰拱形。
如式（10）、式（11）所示，可以看出一個階數對應唯一一個奇異指數和奇異維數，具有唯一性。正常人的奇異譜曲線為單峰拱形，通過觀察采集數據的奇異譜曲線效果對數據進行了篩選，不正常奇異譜曲線如圖 8 所示。數據多重分形奇異譜曲線若出現類似圖 8 的曲線，則刪去該組數據。最終選取 22 組有/無壓力腦電信號數據進行分析。同時對選取的數據進行了統計分析，如表 1、表 2 所示。腦電信號多重分形譜寬度分析結果如表 3 所示。
圖8
不正常奇異譜曲線
Figure8.
Abnormal singular spectrum curve diagram
如表 1 所示，是兩組數據的統計特征不同。如表 2 所示，分析了 Levene 方差齊性檢驗，檢驗兩組方差是否齊，F=0.55，P=0.46 說明方差是齊的。均值方程t 檢驗中P=0.004，說明結果具有統計學意義。
表1
腦電信號統計分析
Table1.
Statistics of two groups of EEG
表2
獨立樣本t 檢驗
Table2.
Independent samplet test
表3
腦電信號多重分形奇異指數
Table3.
Singular index of multi-fractal de-trended fuctuation analysis with electroencephalogram
不同的人對壓力反應不同，所表現出來的多重分形譜不同，即存在個體差異。分析過程中，減去偏差較大的數據，選取 22 組數據進行分析。分別求有/無壓力腦電信號數據奇異譜寬度的最大值、最小值以及平均值并進行比較，比較結果如表 4 所示。
表4
奇異譜寬度比較
Table4.
Comparison of singular spectrum width
有壓力腦電信號的奇異譜寬度最大值、最小值和平均值分別為 0.81、0.35 和 0.59，均高于無壓力腦電信號的奇異譜寬度。基于以上分析，在有效分析得到最優分形階數的基礎上，基于多重分形譜寬度等參數能夠有效分析有/無壓力腦電信號數據，進而有效評估有/無心理壓力狀態。
4 結論
本文基于多重分形去趨勢分析方法實現腦電信號分析與壓力狀態評估。分析了多重分形去趨勢波動分析方法中的質量指數、奇異維數、奇異指數、Hurst 指數與階數的關系。分析結果表明：質量指數和奇異維數隨階數變化不明顯，而 Hurst 指數和奇異指數隨階數變化明顯，心理壓力腦電信號分形階數選擇為 [—5, 5] 時，能夠充分挖掘腦電信號特征。
有壓力腦電信號的奇異譜寬度的最大值和最小值均明顯高于無壓力腦電信號的奇異譜寬度。奇異譜寬度平均值結果顯示，有壓力狀態下，奇異譜寬度平均值明顯高于無壓力狀態下奇異譜寬度平均值，基于奇異譜寬度分析，可以確定壓力狀態。基于多重分形去趨勢波動分析方法，可以有效地識別心理壓力狀態，進而通過相應的干預和預防能夠幫助改善心理壓力狀態，并創建健康的生活方式。
