最大長序列（m 序列）是在高刺激率條件下研究聽覺誘發電位（AEP）線性/非線性成分的一種有效手段。但在實用中如何評估 m 序列的階數選擇對噪聲抑制能力的影響是有待明確的問題。本研究主要采用兩種類型實測數據評估 m 序列噪聲抑制。根據 m 序列互相關法在去卷積計算過程中的兩個計算步驟，即疊加平均和互相關，分別估算不同階數（5～12 階）條件下自發腦電噪聲抑制比，并選擇 7 階和 9 階 m 序列進行非線性 AEP 實驗。結果表明，m 序列對自發腦電的噪聲抑制比與隨機噪聲條件下的理論值完全符合，證實噪聲抑制比主要取決于序列的總長度或總的記錄時間。對于 7 階和 9 階 m 序列的線性/非線性成分進行相似性比較，提示 AEP 成分對階數不敏感。本研究提供了一種選擇 m 序列的更全面的解決方案，可更好地促進基于 m 序列的非線性 AEP 的應用。
引用本文： 陳韻兒, 詹長安, 彭賢, 符秋養, 王濤. 基于最大長序列的聽覺誘發電位噪聲抑制能力分析. 生物醫學工程學雜志, 2018, 35(2): 266-272. doi: 10.7507/1001-5515.201703065 復制
引言
利用最大長序列（maximum length sequence，m 序列）作為聽覺誘發電位（auditory evoked potential，AEP）的一種刺激方案，可同時提取 AEP 中的線性和非線性成分。m 序列是一種偽隨機二值序列，所包括的數值單元個數由序列的階數確定。研究人員通常采用低階 m 序列，利用線性去卷積技術提取快速刺激條件下的線性 AEP[1-2]，或者采用高階 m 序列提取非線性 AEP[3-4]。該過程實質上是一種線性或非線性的去卷積處理，包括兩個重要步驟：首先按照刺激序列的重復次數做疊加平均處理，然后利用互相關法恢復 AEP 的線性或非線性成分[2-4]。按照上述計算過程，m 序列的階數對疊加平均和互相關的噪聲抑制都有影響，在記錄時間不變的情況下這兩種影響是互為補償的。因此，有必要根據處理過程對噪聲的影響和實際測試結果給出 m 序列的選擇依據。
目前 m 序列的噪聲特性研究已有一些可借鑒的工作，由于恢復誘發反應采用的具體算法有所不同，這些研究所得到的具體形式也有差別。比如 Marsh 利用 sub-average 方法來實現 AEP 的提取，并證明其信噪比低于常規方式的誘發反應。若要獲取與常規方式相當的信噪比，需增加平均次數[5]。m 序列用于耳聲發射問題時，Thornton[6-7] 提出一種簡單的改善信噪比的估算公式，和常規方法相比，m 序列可提升 3 dB 的信噪比。Bohorquez 等[8]則在頻域上進行去線性卷積計算，利用 parseval 原理提出了一種 m 序列信噪比的頻域估計方法。
本文根據疊加平均和互相關技術得到一種 m 序列對隨機噪聲抑制比的計算方法，定量說明平均過程和互相關過程在抑制噪聲方面的互補性，并通過自發腦電（electroencephalogram，EEG）和誘發反應的提取實驗，研究階數和數據記錄時間對噪聲抑制以及誘發反應恢復的影響。
1 方法
1.1 基于 m 序列的非線性系統辨識
數學上對于非線性系統的描述一般用 Volterra 級數形式[4]，該模型描述了有限記憶長度非線性系統的 Volterra 核。去卷積計算時，Volterra 核難以估計，可通過正交化過程[9]轉化為具有正交特性的維納核。如果采用高斯白噪聲輸入，Volterra 核或維納核可估計[10]。但對瞬態聲誘發反應來說，更適合采用脈沖型序列作為刺激。Sutter 等[11]利用 m 序列的相關函數為脈沖函數這一性質，提出一種基于互相關的計算方法。Shi 等[3]則進一步將{1，–1}的二元 m 序列拓展到{1，0}形式的 m 脈沖序列，使得誘發刺激序列的生成更加直接便利。
m 序列產生是根據不同階數的本原多項式所確定的移位寄存器實現[12-13]。在給定階數 r 的情況下，本原多項式的個數有限，為不超過 L = 2r – 1 的所有正整數中質數的個數。m 序列包含單元的個數或序列長度為 L。若選定一個本原多項式，可生成一個原型 m 序列。根據實際問題要求，確定刺激率或最小刺激間隔，根據刺激系統時間離散化的要求，確定刺激序列的采樣率，在原型 m 序列單元之間插零稀疏化，調整最終的刺激間隔。若實際刺激采用的 m 序列，最小刺激間隔為 q，則其他刺激間隔為 q 的整數倍，一個周期的 m 序列的實際長度即為 N = L*q。去卷積過程是使用 m 序列對應的一個{+ 1，– 1}的恢復序列進行互相關運算實現[3]。
整個刺激和信號提取過程如圖 1 所示。由于誘發反應信號十分微弱，實驗時 m 脈沖刺激序列會重復刺激聽覺系統，其中序列中“1”表示刺激聲出現時刻。圖 1a 給出的是一個 5 階 m 脈沖序列，所記錄的 EEG 信號如圖 1b 所示。按 m 序列的一個刺激周期（又稱掃程）的起始點為基準疊加平均，得到的平均信號實際上是刺激掃程誘發反應（見圖 1c）。這個平均過程顯然可以抑制非鎖時的信號和隨機噪聲。互相關法提取誘發反應的過程是利用恢復序列 sr（t）和掃程誘發反應 y（t）做循環互相關（見圖 1c—d），得到互相關信號（見圖 1e），即：
 |
該過程是恢復系統成分的一個關鍵步驟，因為所有的脈沖核切片都分布于互相關信號上，核切片的起始位置由位移函數決定，而位移函數是根據所需的核切片個數以及 m 序列的循環右移，通過相等性比較得到[3]。
以上過程中疊加平均和互相關運算分別對隨機噪聲產生抑制作用。假設背景噪聲不存在時間相關性，可分別估計疊加平均和互相關運算兩過程對隨機噪聲的抑制程度[14]。已知平均處理對隨機噪聲幅度的抑制能力為
，其中，K 為平均次數。設噪聲的均方根（root mean square，RMS）值為 σn，原始記錄掃程反應 y（t）中包含 K 個掃程的 EEG 參與疊加平均，則平均過程對噪聲的抑制比為：
 |
其中，σa 為疊加平均后噪聲的 RMS 值，單位為 dB。假定 EEG 實驗記錄時間相同，不同階數 m 序列的稀疏化系數 q 一樣，那么疊加平均次數與 m 序列階數 r 的關系可以表示為：
 |
相對高階 m 序列誘發信號而言，低階 m 序列誘發信號在疊加平均過程中，疊加平均次數多，抑制能力強。疊加平均后的互相關運算將進一步抑制加性噪聲，使得其 RMS 值從 σa 變為 σ?。式（1）的互相關運算是疊加平均過程得到的掃程誘發反應 y（t）與 sr（t）進行卷積實現的。sr（t）包含 L 個交替極性的脈沖（見圖 1d），該卷積過程相當于將 y（t）進行 L 次移位循環疊加。因此卷積后的噪聲 RMS 值將會增至原來的
倍，考慮到式（1）的系數 2/（L + 1），則互相關運算過程對噪聲的抑制程度定義為相關噪聲抑制比：
 |
其中，近似成立條件是 L >> 1，
單位為 dB。綜合式（2）和（4），EEG 整體的噪聲抑制比為：
 |
其中，當給定 q 后，
可以認為是一個常數。式（5）表明 EEG 信號經過疊加平均與互相關運算兩個過程后，整體噪聲水平取決于信號的總長度 N，而 m 序列的階數對整體噪聲抑制比不產生影響。
圖1
一個 5 階 m 序列去卷積過程示意圖
a. m 脈沖序列；b. 腦電信號；c. 平均信號；d. 恢復序列；e. 互相關信號
Figure1. Diagram of the deconvolution process using an m-sequence of order 5a. m-pulse sequence; b. electroencephalogram; c. averaged signal; d. recovery sequence; e. cross-correlation signal
1.2 基于背景腦電的 m 序列仿真實驗
為驗證上述噪聲抑制比計算公式，我們選取階數從 5 到 12 階的 8 個原型 m 序列。設定時間離散的采樣率為 20 kHz，稀疏化系數 q = 40，這樣得到實際刺激序列的最小刺激間隔為 2 ms，最大刺激率為 500 Hz。一個掃程序列的時間長度分別為{62，126，254，510，1 022，2 046，4 094，8 190} ms。假設每個 m 序列實際記錄時間相等，約為 13.5 min。由于序列總長度必須為刺激掃程的整數倍，每個序列的實際長度經調整后略有差別。根據式（3），我們可以得到 5 到 12 階 m 序列的個數分別為 K = {12 800，6 400，3 200，1 600，800，400，200，100}。EEG 信號選擇在 m 序列刺激條件下記錄的自發腦電部分，即刺激停止階段的背景腦電信號。由于這些信號并不包含任何誘發反應的成分，經過相關濾波后信號的強度即可反映 m 序列對噪聲的抑制情況。
1.3 非線性誘發電位實驗
1.3.1 實驗刺激序列的選取
利用 m 序列作為刺激序列時，系統成分在其互相關信號上的分布位置與 m 序列階數有關。低階 m 序列對應的互相關函數的長度相對較短，非線性成分之間會發生重疊，導致失真。刺激序列階數太高時，對實驗設計、數據處理計算會帶來一定的負擔。綜合這兩個情況，我們選取 7 階和 9 階 m 序列作為刺激序列。
1.3.2 實驗方案
受試者為 9 名聽力正常的青年（3 名女性），年齡（21.5 ± 0.85）歲，均為右利手。所有受試者自述無神經系統、聽覺系統疾病。本文實驗得到南方醫科大學附屬南方醫院倫理委員會的批準，并且各受試者均對實驗知情，同意進行實驗。EEG 采集在聲電屏蔽室進行，采用澳大利亞 Neuroscan 公司的 SynAmps2 系統。實驗前皮膚處理：75% 的醫用酒精拭擦局部皮膚，之后用磨砂膏去除局部角質，再用相同濃度的醫用酒精消毒去污。記錄電極放置于額頭正中發際線處，參考電極置于右耳乳突，眉心接地，并使得阻抗 < 5 kΩ。記錄數據時，受試者放松平躺在椅子上，盡量避免身體移動。記錄的參數設置為：采樣頻率為 20 kHz，帶通濾波范圍為 100～2 000 Hz，測試采用 7 階和 9 階 m 序列，稀疏化系數、采樣頻率和最小刺激間隔的設定與仿真實驗一致。刺激聲為 0.2 ms 持續時間的疏波短聲（click），峰值聲強校準為 82 dB SPL，通過 Etymotic Research ER-3A 耳機進行播放。7 階和 9 階 m 序列記錄的掃程個數分別為 2 000 和 500，對應的實驗時間分別為 8.47 min 和 8.52 min。
2 結果
2.1 自發背景腦電仿真結果
圖 2a 中給出記錄時間約為 13.44 min 時，15 次仿真計算的噪聲抑制比的平均值和標準差（見圖 2a 實線），這些實驗計算值與式（2）給出隨機噪聲條件下的理論值（見圖 2a 虛線）高度一致，幾近重合。對平均信號的互相關運算結果如圖 2b 所示，其中實驗計算值與 m 序列的階數 r 呈線性遞減關系（見圖 2b 實線），與公式（4）的理論值（見圖 2b 虛線）也高度一致。
圖2
噪聲抑制比（ηa 和 η?）與階數 r 的線性關系
a.
a.
圖 3 為選擇 4 種記錄時間條件下，整個過程的噪聲抑制比和 m 序列階數 r 的關系，計算結果與公式（5）給出的理論值一致，并且對于不同記錄時間，整個過程的總體噪聲抑制比是上述疊加平均和互相關運算兩過程的噪聲抑制比之和，與階數無關。這說明不同階數 m 序列的疊加平均和互相關運算兩過程在噪聲抑制方面具有互補性。記錄時間相同時，兩者之和在不同階數情況下近似不變。因為在此情況下，單次的平均運算和互相關過程中加法運算對隨機噪聲具有相同的抑制能力，而兩者的運算次數存在此消彼長的關系。由于總的記錄時間和一個掃程的長度存在倍數關系，實際上每個階數對應的記錄時間和設定值略有差別，這可能是導致圖 3 中實驗值（實線）偏離理論值（虛線）較大的原因。
圖3
不同記錄時間條件下，整個去卷積過程的噪聲抑制比 ηa? 與階數 r 的關系
Figure3.
Relationships between the whole noise attenuation ration ηa? and m-sequence order r for different recording times
2.2 m 序列實測實驗結果
在 AEP 實測實驗時，EEG 信號除了各種偽跡和噪聲外還包括相對固定的確定性誘發反應，即 AEP 成分。經過疊加平均和互相關運算兩過程后，噪聲能量被抑制到一定程度時，AEP 成分對噪聲抑制比計算的影響會導致計算誤差。一種簡便的消除誘發反應信號影響的方法是奇偶交替平均技術[15]，即將所有參與平均的 EEG 樣本的奇數平均或偶數平均取反后再平均，這樣可抵消其中鎖相的誘發反應。圖 4 分別給出 7 階和 9 階所有受試者不同次數平均處理的噪聲抑制比的結果。可見直接計算時，誘發反應的影響使得噪聲抑制比（‘×’）普遍大于理論值（虛線），且波動較大。利用奇偶交替平均技術消除誘發反應影響后，所有受試者的結果（‘●’）與理論值都很好地吻合了。本例中 7 階和 9 階最多的疊加次數分別為 2 000 和 500，其對應時間大約為 8.5 min。
圖4
疊加平均過程的噪聲抑制比與疊加次數的關系
Figure4.
Relationship between the noise attenuation ratio of averaging process and sweep number for 7 and 9 m-sequence orders
如果線性和非線性 AEP 成分存在，它們將分布在互相關信號的不同部位，并或強或弱地突出于背景噪聲。除了利用奇偶交替技術抵消 AEP 成分外，還可利用位移函數定位并剪除這些有效 AEP 成分，僅保留殘余噪聲部分。圖 5a-c 反映了疊加平均與互相關運算兩個過程的總體噪聲抑制比。其中水平虛線是式（5）給出兩個過程整體的理論值，因為實驗時間近似，所以 7 階和 9 階 m 序列誘發的 EEG 總體噪聲抑制比的理論值近似相同。如果不排除有效的 AEP 成分，采用全部信息計算整體噪聲抑制比，不管是 7 階（‘-×’）還是 9 階（‘-○’）對應的結果均與理論值偏離較大（見圖 5a）；將互相關信號上的 AEP 成分剪除后，7 階和 9 階對應的總體噪聲抑制比結果與理論值吻合（見圖 5b）。類似地利用奇偶交替技術抵消 AEP 成分，7 階和 9 階對應的總體噪聲抑制比與理論值也基本吻合（見圖 5c）。這些結果顯示，平均與互相關處理去除確定性誘發反應后，噪聲抑制比與理論值吻合較好。
圖5
兩個 m 序列在相同實驗時間條件下，每個受試者的總體噪聲抑制比
a. 含有誘發反應；b. 去除互相關階段誘發反應成分；c. 交替平均去除誘發反應成分；虛線表示式（5）計算得到的理論值
Figure5. Individual noise attenuation ratios of two m-sequences for true data with different noise approximation methodsa. noise with evoked responses; b. response excluded from the correlation process; c. noise-only by alternative averaging; dashed lines indicate theoretical values given by equation (5)
2.3 兩個 m 序列對應的線性和非線性成分的相似性比較
本節旨在比較在記錄時間近似的情況下，兩個階數的 m 序列所誘發的 AEP 成分是否相近，從而初步判斷階數是否對有效 AEP 成分產生影響。圖 6 給出 7 階（細實線）和 9 階（粗實線）m 序列得到的 4 個 AEP 的波形。其中包括一個線性 AEP 成分，即一階核切片 KS11，以及三個非線性 AEP 成分：二階核切片 KS21、KS22 和三階核切片 KS31。選擇系統的記憶長度為 T = 20 ms。視覺評估可見，兩個 m 序列的誘發反應相似性較好，特別是 KS11 中標出的三個主波對應的幅度和潛伏期一致性很高。而在當前信噪比條件下，三個非線性成分的三個主波難以辨識。
圖6
兩個 m 序列誘發的線性和非線性成分對比
Figure6.
Waveform comparison of the AEP components for two m-sequences
圖 7 是每個受試者的 AEP 成分的總體平均結果。圖中 KS11 顯示出聽覺腦干反應（auditory brainstem response，ABR）的幾個主要成分，如圖中所標的Ⅰ、Ⅲ、Ⅴ波。由于 m 序列最小間隔為 2 ms，即最大刺激率為 500 Hz，這些成分的潛伏期與常規記錄方式相比數值較大，與一般觀測結果相符[1, 16]，說明兩個刺激序列可成功獲得線性 ABR 成分，并且具有良好的一致性。對于其余三個非線性 AEP 成分，不同刺激序列誘發的 AEP 成分波形之間能夠較好地吻合。為了量化兩個刺激序列誘發的 AEP 成分波形之間的相似度，分別計算相同 AEP 成分波形之間的相關系數，其中線性 AEP 成分對應為 0.95，三個非線性 AEP 成分分別對應為{0.81，0.53，0.79}。
圖7
兩個 m 序列誘發的 AEP 各成分總體平均波形比較
Figure7.
Comparison of the grand averaged AEP components for two m-sequences
3 討論
對于常規方法而言，信噪比的提高與疊加刺激平方根呈正相關，如果 m 序列方法想要獲得與常規方法相當的信噪比，需要增加疊加次數，相當于延長實驗時間[5]。但是高頻刺激聲可以增加聽覺系統的負荷，提高聽覺系統和腦部病變檢測的靈敏度，使得檢測結果更為準確[17]。因此利用 m 序列作為高頻刺激序列對于研究 AEP 線性和非線性特性具有一定的研究價值。由于 m 序列階數的選擇會對其誘發的 AEP 線性和非線性成分產生影響，以及 m 序列階數對噪聲抑制能力的評估尚無明確的評判標準，所以本文通過仿真和實測實驗研究了 m 序列階數對噪聲抑制的影響，同時利用實測實驗研究 m 序列階數對 AEP 成分的影響。
本文根據疊加平均以及互相關運算對噪聲抑制程度的計算公式，采用自發腦電和實際 m 序列誘發 AEP 實驗進行驗證。表明 AEP 的噪聲抑制能力僅取決于總的記錄時間，與 m 序列的階數無關。理論上只要記錄時間相近，任意階數 m 序列的噪聲抑制情況都近似，而且給定記錄時間，便可預判任意階數 m 序列的噪聲抑制程度，或者根據噪聲抑制程度可以反推所需的記錄時間。此外，本文還分析了 7 階和 9 階 m 序列誘發的線性/非線性 AEP 成分之間的相似性，表明在刺激率相同時，兩個序列的有效 AEP 成分波形相近，提示階數對 AEP 成分亦無影響。但是，m 序列的階數以及相同階數的不同 m 序列可能在其他方面影響提取誘發反應的性能。比如 AEP 成分的分布位置也與 m 序列的性質有關，可能會導致 AEP 成分的重疊失真[17]，提示在 m 序列的選擇方面還需要全面考慮。
4 結論
實測數據的 m 序列的噪聲抑制比與理論值高度相符。在相同記錄時間條件下，m 序列的階數對誘發反應的噪聲抑制能力沒有影響。結合 m 序列對非線性成分的分布性質影響的研究結果[18]，本研究可提供一種選擇 m 序列的更全面的解決方案。
致謝：感謝南方醫科大學林霖老師對本文公式的推導給予了幫助。
引言
利用最大長序列（maximum length sequence，m 序列）作為聽覺誘發電位（auditory evoked potential，AEP）的一種刺激方案，可同時提取 AEP 中的線性和非線性成分。m 序列是一種偽隨機二值序列，所包括的數值單元個數由序列的階數確定。研究人員通常采用低階 m 序列，利用線性去卷積技術提取快速刺激條件下的線性 AEP[1-2]，或者采用高階 m 序列提取非線性 AEP[3-4]。該過程實質上是一種線性或非線性的去卷積處理，包括兩個重要步驟：首先按照刺激序列的重復次數做疊加平均處理，然后利用互相關法恢復 AEP 的線性或非線性成分[2-4]。按照上述計算過程，m 序列的階數對疊加平均和互相關的噪聲抑制都有影響，在記錄時間不變的情況下這兩種影響是互為補償的。因此，有必要根據處理過程對噪聲的影響和實際測試結果給出 m 序列的選擇依據。
目前 m 序列的噪聲特性研究已有一些可借鑒的工作，由于恢復誘發反應采用的具體算法有所不同，這些研究所得到的具體形式也有差別。比如 Marsh 利用 sub-average 方法來實現 AEP 的提取，并證明其信噪比低于常規方式的誘發反應。若要獲取與常規方式相當的信噪比，需增加平均次數[5]。m 序列用于耳聲發射問題時，Thornton[6-7] 提出一種簡單的改善信噪比的估算公式，和常規方法相比，m 序列可提升 3 dB 的信噪比。Bohorquez 等[8]則在頻域上進行去線性卷積計算，利用 parseval 原理提出了一種 m 序列信噪比的頻域估計方法。
本文根據疊加平均和互相關技術得到一種 m 序列對隨機噪聲抑制比的計算方法，定量說明平均過程和互相關過程在抑制噪聲方面的互補性，并通過自發腦電（electroencephalogram，EEG）和誘發反應的提取實驗，研究階數和數據記錄時間對噪聲抑制以及誘發反應恢復的影響。
1 方法
1.1 基于 m 序列的非線性系統辨識
數學上對于非線性系統的描述一般用 Volterra 級數形式[4]，該模型描述了有限記憶長度非線性系統的 Volterra 核。去卷積計算時，Volterra 核難以估計，可通過正交化過程[9]轉化為具有正交特性的維納核。如果采用高斯白噪聲輸入，Volterra 核或維納核可估計[10]。但對瞬態聲誘發反應來說，更適合采用脈沖型序列作為刺激。Sutter 等[11]利用 m 序列的相關函數為脈沖函數這一性質，提出一種基于互相關的計算方法。Shi 等[3]則進一步將{1，–1}的二元 m 序列拓展到{1，0}形式的 m 脈沖序列，使得誘發刺激序列的生成更加直接便利。
m 序列產生是根據不同階數的本原多項式所確定的移位寄存器實現[12-13]。在給定階數 r 的情況下，本原多項式的個數有限，為不超過 L = 2r – 1 的所有正整數中質數的個數。m 序列包含單元的個數或序列長度為 L。若選定一個本原多項式，可生成一個原型 m 序列。根據實際問題要求，確定刺激率或最小刺激間隔，根據刺激系統時間離散化的要求，確定刺激序列的采樣率，在原型 m 序列單元之間插零稀疏化，調整最終的刺激間隔。若實際刺激采用的 m 序列，最小刺激間隔為 q，則其他刺激間隔為 q 的整數倍，一個周期的 m 序列的實際長度即為 N = L*q。去卷積過程是使用 m 序列對應的一個{+ 1，– 1}的恢復序列進行互相關運算實現[3]。
整個刺激和信號提取過程如圖 1 所示。由于誘發反應信號十分微弱，實驗時 m 脈沖刺激序列會重復刺激聽覺系統，其中序列中“1”表示刺激聲出現時刻。圖 1a 給出的是一個 5 階 m 脈沖序列，所記錄的 EEG 信號如圖 1b 所示。按 m 序列的一個刺激周期（又稱掃程）的起始點為基準疊加平均，得到的平均信號實際上是刺激掃程誘發反應（見圖 1c）。這個平均過程顯然可以抑制非鎖時的信號和隨機噪聲。互相關法提取誘發反應的過程是利用恢復序列 sr（t）和掃程誘發反應 y（t）做循環互相關（見圖 1c—d），得到互相關信號（見圖 1e），即：
|
該過程是恢復系統成分的一個關鍵步驟，因為所有的脈沖核切片都分布于互相關信號上，核切片的起始位置由位移函數決定，而位移函數是根據所需的核切片個數以及 m 序列的循環右移，通過相等性比較得到[3]。
以上過程中疊加平均和互相關運算分別對隨機噪聲產生抑制作用。假設背景噪聲不存在時間相關性，可分別估計疊加平均和互相關運算兩過程對隨機噪聲的抑制程度[14]。已知平均處理對隨機噪聲幅度的抑制能力為
，其中，K 為平均次數。設噪聲的均方根（root mean square，RMS）值為 σn，原始記錄掃程反應 y（t）中包含 K 個掃程的 EEG 參與疊加平均，則平均過程對噪聲的抑制比為：
|
其中，σa 為疊加平均后噪聲的 RMS 值，單位為 dB。假定 EEG 實驗記錄時間相同，不同階數 m 序列的稀疏化系數 q 一樣，那么疊加平均次數與 m 序列階數 r 的關系可以表示為：
|
相對高階 m 序列誘發信號而言，低階 m 序列誘發信號在疊加平均過程中，疊加平均次數多，抑制能力強。疊加平均后的互相關運算將進一步抑制加性噪聲，使得其 RMS 值從 σa 變為 σ?。式（1）的互相關運算是疊加平均過程得到的掃程誘發反應 y（t）與 sr（t）進行卷積實現的。sr（t）包含 L 個交替極性的脈沖（見圖 1d），該卷積過程相當于將 y（t）進行 L 次移位循環疊加。因此卷積后的噪聲 RMS 值將會增至原來的
倍，考慮到式（1）的系數 2/（L + 1），則互相關運算過程對噪聲的抑制程度定義為相關噪聲抑制比：
|
其中，近似成立條件是 L >> 1，
單位為 dB。綜合式（2）和（4），EEG 整體的噪聲抑制比為：
|
其中，當給定 q 后，
可以認為是一個常數。式（5）表明 EEG 信號經過疊加平均與互相關運算兩個過程后，整體噪聲水平取決于信號的總長度 N，而 m 序列的階數對整體噪聲抑制比不產生影響。
圖1
一個 5 階 m 序列去卷積過程示意圖
a. m 脈沖序列；b. 腦電信號；c. 平均信號；d. 恢復序列；e. 互相關信號
Figure1. Diagram of the deconvolution process using an m-sequence of order 5a. m-pulse sequence; b. electroencephalogram; c. averaged signal; d. recovery sequence; e. cross-correlation signal
1.2 基于背景腦電的 m 序列仿真實驗
為驗證上述噪聲抑制比計算公式，我們選取階數從 5 到 12 階的 8 個原型 m 序列。設定時間離散的采樣率為 20 kHz，稀疏化系數 q = 40，這樣得到實際刺激序列的最小刺激間隔為 2 ms，最大刺激率為 500 Hz。一個掃程序列的時間長度分別為{62，126，254，510，1 022，2 046，4 094，8 190} ms。假設每個 m 序列實際記錄時間相等，約為 13.5 min。由于序列總長度必須為刺激掃程的整數倍，每個序列的實際長度經調整后略有差別。根據式（3），我們可以得到 5 到 12 階 m 序列的個數分別為 K = {12 800，6 400，3 200，1 600，800，400，200，100}。EEG 信號選擇在 m 序列刺激條件下記錄的自發腦電部分，即刺激停止階段的背景腦電信號。由于這些信號并不包含任何誘發反應的成分，經過相關濾波后信號的強度即可反映 m 序列對噪聲的抑制情況。
1.3 非線性誘發電位實驗
1.3.1 實驗刺激序列的選取
利用 m 序列作為刺激序列時，系統成分在其互相關信號上的分布位置與 m 序列階數有關。低階 m 序列對應的互相關函數的長度相對較短，非線性成分之間會發生重疊，導致失真。刺激序列階數太高時，對實驗設計、數據處理計算會帶來一定的負擔。綜合這兩個情況，我們選取 7 階和 9 階 m 序列作為刺激序列。
1.3.2 實驗方案
受試者為 9 名聽力正常的青年（3 名女性），年齡（21.5 ± 0.85）歲，均為右利手。所有受試者自述無神經系統、聽覺系統疾病。本文實驗得到南方醫科大學附屬南方醫院倫理委員會的批準，并且各受試者均對實驗知情，同意進行實驗。EEG 采集在聲電屏蔽室進行，采用澳大利亞 Neuroscan 公司的 SynAmps2 系統。實驗前皮膚處理：75% 的醫用酒精拭擦局部皮膚，之后用磨砂膏去除局部角質，再用相同濃度的醫用酒精消毒去污。記錄電極放置于額頭正中發際線處，參考電極置于右耳乳突，眉心接地，并使得阻抗 < 5 kΩ。記錄數據時，受試者放松平躺在椅子上，盡量避免身體移動。記錄的參數設置為：采樣頻率為 20 kHz，帶通濾波范圍為 100～2 000 Hz，測試采用 7 階和 9 階 m 序列，稀疏化系數、采樣頻率和最小刺激間隔的設定與仿真實驗一致。刺激聲為 0.2 ms 持續時間的疏波短聲（click），峰值聲強校準為 82 dB SPL，通過 Etymotic Research ER-3A 耳機進行播放。7 階和 9 階 m 序列記錄的掃程個數分別為 2 000 和 500，對應的實驗時間分別為 8.47 min 和 8.52 min。
2 結果
2.1 自發背景腦電仿真結果
圖 2a 中給出記錄時間約為 13.44 min 時，15 次仿真計算的噪聲抑制比的平均值和標準差（見圖 2a 實線），這些實驗計算值與式（2）給出隨機噪聲條件下的理論值（見圖 2a 虛線）高度一致，幾近重合。對平均信號的互相關運算結果如圖 2b 所示，其中實驗計算值與 m 序列的階數 r 呈線性遞減關系（見圖 2b 實線），與公式（4）的理論值（見圖 2b 虛線）也高度一致。
圖2
噪聲抑制比（ηa 和 η?）與階數 r 的線性關系
a.
a.
圖 3 為選擇 4 種記錄時間條件下，整個過程的噪聲抑制比和 m 序列階數 r 的關系，計算結果與公式（5）給出的理論值一致，并且對于不同記錄時間，整個過程的總體噪聲抑制比是上述疊加平均和互相關運算兩過程的噪聲抑制比之和，與階數無關。這說明不同階數 m 序列的疊加平均和互相關運算兩過程在噪聲抑制方面具有互補性。記錄時間相同時，兩者之和在不同階數情況下近似不變。因為在此情況下，單次的平均運算和互相關過程中加法運算對隨機噪聲具有相同的抑制能力，而兩者的運算次數存在此消彼長的關系。由于總的記錄時間和一個掃程的長度存在倍數關系，實際上每個階數對應的記錄時間和設定值略有差別，這可能是導致圖 3 中實驗值（實線）偏離理論值（虛線）較大的原因。
圖3
不同記錄時間條件下，整個去卷積過程的噪聲抑制比 ηa? 與階數 r 的關系
Figure3.
Relationships between the whole noise attenuation ration ηa? and m-sequence order r for different recording times
2.2 m 序列實測實驗結果
在 AEP 實測實驗時，EEG 信號除了各種偽跡和噪聲外還包括相對固定的確定性誘發反應，即 AEP 成分。經過疊加平均和互相關運算兩過程后，噪聲能量被抑制到一定程度時，AEP 成分對噪聲抑制比計算的影響會導致計算誤差。一種簡便的消除誘發反應信號影響的方法是奇偶交替平均技術[15]，即將所有參與平均的 EEG 樣本的奇數平均或偶數平均取反后再平均，這樣可抵消其中鎖相的誘發反應。圖 4 分別給出 7 階和 9 階所有受試者不同次數平均處理的噪聲抑制比的結果。可見直接計算時，誘發反應的影響使得噪聲抑制比（‘×’）普遍大于理論值（虛線），且波動較大。利用奇偶交替平均技術消除誘發反應影響后，所有受試者的結果（‘●’）與理論值都很好地吻合了。本例中 7 階和 9 階最多的疊加次數分別為 2 000 和 500，其對應時間大約為 8.5 min。
圖4
疊加平均過程的噪聲抑制比與疊加次數的關系
Figure4.
Relationship between the noise attenuation ratio of averaging process and sweep number for 7 and 9 m-sequence orders
如果線性和非線性 AEP 成分存在，它們將分布在互相關信號的不同部位，并或強或弱地突出于背景噪聲。除了利用奇偶交替技術抵消 AEP 成分外，還可利用位移函數定位并剪除這些有效 AEP 成分，僅保留殘余噪聲部分。圖 5a-c 反映了疊加平均與互相關運算兩個過程的總體噪聲抑制比。其中水平虛線是式（5）給出兩個過程整體的理論值，因為實驗時間近似，所以 7 階和 9 階 m 序列誘發的 EEG 總體噪聲抑制比的理論值近似相同。如果不排除有效的 AEP 成分，采用全部信息計算整體噪聲抑制比，不管是 7 階（‘-×’）還是 9 階（‘-○’）對應的結果均與理論值偏離較大（見圖 5a）；將互相關信號上的 AEP 成分剪除后，7 階和 9 階對應的總體噪聲抑制比結果與理論值吻合（見圖 5b）。類似地利用奇偶交替技術抵消 AEP 成分，7 階和 9 階對應的總體噪聲抑制比與理論值也基本吻合（見圖 5c）。這些結果顯示，平均與互相關處理去除確定性誘發反應后，噪聲抑制比與理論值吻合較好。
圖5
兩個 m 序列在相同實驗時間條件下，每個受試者的總體噪聲抑制比
a. 含有誘發反應；b. 去除互相關階段誘發反應成分；c. 交替平均去除誘發反應成分；虛線表示式（5）計算得到的理論值
Figure5. Individual noise attenuation ratios of two m-sequences for true data with different noise approximation methodsa. noise with evoked responses; b. response excluded from the correlation process; c. noise-only by alternative averaging; dashed lines indicate theoretical values given by equation (5)
2.3 兩個 m 序列對應的線性和非線性成分的相似性比較
本節旨在比較在記錄時間近似的情況下，兩個階數的 m 序列所誘發的 AEP 成分是否相近，從而初步判斷階數是否對有效 AEP 成分產生影響。圖 6 給出 7 階（細實線）和 9 階（粗實線）m 序列得到的 4 個 AEP 的波形。其中包括一個線性 AEP 成分，即一階核切片 KS11，以及三個非線性 AEP 成分：二階核切片 KS21、KS22 和三階核切片 KS31。選擇系統的記憶長度為 T = 20 ms。視覺評估可見，兩個 m 序列的誘發反應相似性較好，特別是 KS11 中標出的三個主波對應的幅度和潛伏期一致性很高。而在當前信噪比條件下，三個非線性成分的三個主波難以辨識。
圖6
兩個 m 序列誘發的線性和非線性成分對比
Figure6.
Waveform comparison of the AEP components for two m-sequences
圖 7 是每個受試者的 AEP 成分的總體平均結果。圖中 KS11 顯示出聽覺腦干反應（auditory brainstem response，ABR）的幾個主要成分，如圖中所標的Ⅰ、Ⅲ、Ⅴ波。由于 m 序列最小間隔為 2 ms，即最大刺激率為 500 Hz，這些成分的潛伏期與常規記錄方式相比數值較大，與一般觀測結果相符[1, 16]，說明兩個刺激序列可成功獲得線性 ABR 成分，并且具有良好的一致性。對于其余三個非線性 AEP 成分，不同刺激序列誘發的 AEP 成分波形之間能夠較好地吻合。為了量化兩個刺激序列誘發的 AEP 成分波形之間的相似度，分別計算相同 AEP 成分波形之間的相關系數，其中線性 AEP 成分對應為 0.95，三個非線性 AEP 成分分別對應為{0.81，0.53，0.79}。
圖7
兩個 m 序列誘發的 AEP 各成分總體平均波形比較
Figure7.
Comparison of the grand averaged AEP components for two m-sequences
3 討論
對于常規方法而言，信噪比的提高與疊加刺激平方根呈正相關，如果 m 序列方法想要獲得與常規方法相當的信噪比，需要增加疊加次數，相當于延長實驗時間[5]。但是高頻刺激聲可以增加聽覺系統的負荷，提高聽覺系統和腦部病變檢測的靈敏度，使得檢測結果更為準確[17]。因此利用 m 序列作為高頻刺激序列對于研究 AEP 線性和非線性特性具有一定的研究價值。由于 m 序列階數的選擇會對其誘發的 AEP 線性和非線性成分產生影響，以及 m 序列階數對噪聲抑制能力的評估尚無明確的評判標準，所以本文通過仿真和實測實驗研究了 m 序列階數對噪聲抑制的影響，同時利用實測實驗研究 m 序列階數對 AEP 成分的影響。
本文根據疊加平均以及互相關運算對噪聲抑制程度的計算公式，采用自發腦電和實際 m 序列誘發 AEP 實驗進行驗證。表明 AEP 的噪聲抑制能力僅取決于總的記錄時間，與 m 序列的階數無關。理論上只要記錄時間相近，任意階數 m 序列的噪聲抑制情況都近似，而且給定記錄時間，便可預判任意階數 m 序列的噪聲抑制程度，或者根據噪聲抑制程度可以反推所需的記錄時間。此外，本文還分析了 7 階和 9 階 m 序列誘發的線性/非線性 AEP 成分之間的相似性，表明在刺激率相同時，兩個序列的有效 AEP 成分波形相近，提示階數對 AEP 成分亦無影響。但是，m 序列的階數以及相同階數的不同 m 序列可能在其他方面影響提取誘發反應的性能。比如 AEP 成分的分布位置也與 m 序列的性質有關，可能會導致 AEP 成分的重疊失真[17]，提示在 m 序列的選擇方面還需要全面考慮。
4 結論
實測數據的 m 序列的噪聲抑制比與理論值高度相符。在相同記錄時間條件下，m 序列的階數對誘發反應的噪聲抑制能力沒有影響。結合 m 序列對非線性成分的分布性質影響的研究結果[18]，本研究可提供一種選擇 m 序列的更全面的解決方案。
致謝：感謝南方醫科大學林霖老師對本文公式的推導給予了幫助。
