本文基于已有的卷積核補償(CKC)方法, 提出了一種新的信號分解方法。該方法與自組織映射神經網絡相結合, 首先找出一個在某一時刻具有發放活動的脈沖序列, 其次對這個脈沖序列的一些較大值所對應的時刻利用自組織映射神經網絡進行分類, 然后利用分類后的時刻所對應的測量信號的值求出最終的一個信號源的發放序列。通過隨機混合矩陣合成產生的仿真信號進行測試, 表明所提出的方法是有效的。
引用本文： 寧勇, 朱善安, 趙玉明. 基于自組織映射神經網絡和卷積核補償的多通道表面肌電信號的盲源分離方法. 生物醫學工程學雜志, 2015, 32(1): 1-7. doi: 10.7507/1001-5515.20150001 復制
引言
生物醫學信號往往包含著復雜的信息，通常由多個信號源的信號綜合所組成，比如神經元和神經纖維。這些信號源通常被看作是統計獨立的，或者至少是弱相關的，但是它們的混合過程實際上是不知道的，因此需要對所獲得的這類信號進行盲分解。
過去一段時間以來針對肌電(electromyogram, EMG)信號的研究是很活躍的，因為無論其在肌肉康復治療中還是在臨床診斷中都有著重要的現實意義，研究者們還通過測量EMG信號研究神經肌肉的控制機理以及在其它方面進行應用，例如假肢控制[1-2]，手語識別[3]等。EMG信號中的動作電位(action potential, AP)通常由多個運動單位(motor unit, MU)所產生的，其被看作是隨機干擾的非平穩信號，具有復雜的特性，尤其是表面肌電(surface electromyogram, sEMG)信號的復雜程度更要高一些，不但其測量信號的電極與所接觸皮膚的接觸特性不良好，而且還有各種噪聲的干擾，因此加劇了對sEMG信號分解和提取臨床相關信息的困難。
近年來，高密度陣列式的sEMG電極被用來進行sEMG信號的檢測，針對這種檢測方式，一些sEMG信號處理的方法被陸續提出，包括模式識別的方法和盲源分離的方法。Kleine等[4]指出二維空間濾波器在sEMG信號分解中對檢測發放精度的重要性。Gazzoni等[5]采用分類技術和模板匹配方法對信號進行分解研究。Wood等[6]通過有限元分析方法來研究運動單元的深度和皮膚表面電位分布之間的關系。Garcia等[7]把sEMG信號的產生看作是運動單元動作電位(motor unit action potential, MUAP)序列的瞬時混合過程。Gligorijevic等[8]先對運動單元電位波形進行聚類，然后再用一種分解技術對信號進行分解。Holobar等[9-13]提出了卷積核補償(gradient convolution kernel compensation, CKC)方法和梯度優化的CKC方法來處理卷積混合的運動單元發放源信號，這種方法不需要計算未知的混合矩陣并且對發放源信號直接進行計算。CKC方法被證明是很有效的，如果sEMG信號的質量較好，可以同時重建出30多個發放源[10]。
本文在現有CKC方法的理論研究和實踐應用的基礎之上，嘗試提出一種新的方法，該方法將結合自組織映射神經網絡對信號進行處理，以期對現有的CKC方法的計算過程進行進一步的簡化并獲得更好的處理結果。
1 原理與方法
1.1 數據模型和卷積核補償
對于一個線性時不變(linear time invariable, LTI)的多輸入多輸出(multi-input multi-output, MIMO)系統，CKC的數學模型可以用以下方程表示[9-10]：
| $ \boldsymbol{X}\text{ }(n)=\boldsymbol{H}\bar{s}\left(n \right)+\boldsymbol{e}\left(n \right), $ |
其中X(n)=[x1(n)，…, xM(n)]T為卷積后的多通道檢測信號向量，xj(n)為第j個通道中的第n個采樣值，e(n)=[e1(n), …, eM(n)]T為零均值白噪聲向量，
| $ \boldsymbol{H}=\left(\begin{align} & {{h}_{11}}\left(0 \right)\cdots {{h}_{11}}\left(P-1 \right)\ \ \ {{h}_{12}}\left(0 \right)\cdots {{h}_{12}}\left(P-1 \right)\cdots {{h}_{1N}}\left(0 \right)\cdots {{h}_{1N}}\left(P-1 \right)\\ & {{h}_{21}}\left(0 \right)\cdots {{h}_{21}}\left(P-1 \right)\ \ \ {{h}_{22}}\left(0 \right)\cdots {{h}_{22}}\left(P-1 \right)\cdots {{h}_{2N}}\left(0 \right)\cdots {{h}_{2N}}\left(P-1 \right)\\ & \ \ \ \vdots \ \ \ \ \ \cdots \ \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \ \ \cdots \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \cdots \ \ \ \ \vdots \ \ \ \ \ \ \ \cdots \ \ \ \ \ \ \ \ \ \ \vdots \\ & {{h}_{M1}}\left(0 \right)\cdots {{h}_{M1}}\left(P-1 \right)\ \ \ {{h}_{M2}}\left(0 \right)\cdots {{h}_{M2}}\left(P-1 \right)\cdots {{h}_{MN}}\left(0 \right)\cdots {{h}_{MN}}\left(P-1 \right)\\ \end{align} \right), $ |
其中hij是在第i個通道中第j個信號源的源脈沖序列的響應，si(n)為第i個信號源的發放序列，可以用以下脈沖序列的模型表示：
| $ {{s}_{i}}\left(n \right)=\sum\limits_{-\infty }^{\infty }{\delta \left[n-{{S}_{i}}\left(k \right)\right]}, $ |
其中δ(·)為狄拉克沖激函數，Si(k)為第i個信號源的第k個發放時刻。
為了改善通道數和未知信號源的數目之間的比例，以提高分解效果，在模型(1)中，對每一個X(n)分量進行k-1個延遲擴展，即
| $ \begin{align} & \gamma \left(n \right)={{{\bar{\boldsymbol{X}}}}^{\boldsymbol{T}}}\left(n \right)\boldsymbol{C}_{\bar{\boldsymbol{X}}\bar{\boldsymbol{X}}}^{-\ 1}\bar{\boldsymbol{X}}\left(n \right)=\bar{s}\left(n \right){{\boldsymbol{H}}^{T}}{{\left({{\boldsymbol{H}}^{T}} \right)}^{-1}}\boldsymbol{C}_{\bar{s}\bar{s}}^{-1}{{\boldsymbol{H}}^{-1}} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \boldsymbol{H}\ \bar{s}\left(n \right)={{{\bar{s}}}^{T}}\left(n \right)\boldsymbol{C}_{\bar{s}\bar{s}}^{-1}\bar{s}\left(n \right), \\ \end{align} $ |
其中C
如果式(4)中右側第一項用
| $ {{s}_{{{n}_{0}}}}\left(n \right)={{{\bar{\boldsymbol{X}}}}^{T}}\left({{n}_{0}} \right)\boldsymbol{C}_{\bar{\boldsymbol{X}}\bar{\boldsymbol{X}}}^{-\ 1}\bar{\boldsymbol{X}}\left(n \right), $ |
已經證明如果第j個信號源在時刻n0發放，sn0(n)可表示為對第j個信號源發放序列的估計[9]，
| $ {{s}_{{{n}_{0}}}}\left(n \right)\approx {{{\bar{s}}}_{j}}\left(n \right), $ |
這樣可以通過式(5)計算出一個信號源的脈沖發放序列，以進行信號的分解。為了改善sn0(n)的計算質量，式(5)中的
| $ {{\boldsymbol{C}}_{{{\boldsymbol{X}}_{sj}}}}=\frac{1}{card\left({{\varphi }_{n}} \right)}\sum{\bar{\boldsymbol{X}}\left({{\varphi }_{n}} \right)}, $ |
其中φn={n0, n1, …, nm}，card(φn)表示φn中所包含元素的數目。φn中所有的時刻應該被同一個信號源的脈沖序列的發放時刻所包括，或是絕大部分被包括，這樣計算出的效果會更好。
1.2 基于自組織映射神經網絡的卷積核補償方法
自組織映射神經網絡[14-16]是基于無監督學習方法的神經網絡的一種重要類型。自組織映射網絡理論最早是由芬蘭赫爾辛基理工大學的Kohohen于1981年提出的。此后，伴隨著神經網絡在20世紀80年代中后期的迅速發展，自組織映射理論及其應用也有了長足的進步。
自組織影射神經網絡的訓練過程是無監督的聚類過程, 其通過對輸入樣本的反復學習, 使連接權向量的空間分布與輸入樣本的概率分布趨于一致, 進而將多維的任意輸入模式在輸出層映射成一維、二維或多維離散圖形, 并且圖形保持輸入的拓撲結構不變。自組織映射網絡的學習算法抽象意義上包括：競爭、合作和更新三個階段。
(1)競爭：對于每一個輸入樣本, 網絡中的神經元對它們各自的判別函數值進行計算, 這個判別函數值對神經元之間的競爭提供了基礎, 最終具有函數最大值的神經元將成為競爭的勝利者。
(2)合作：獲勝的神經元決定了興奮神經元的拓撲鄰域空間位置, 進而為鄰域神經元提供合作的基礎。
(3)更新：興奮神經元通過對它們各自權值的適當調節以增強該輸入模式的判別函數值, 使以后相似模式輸入的響應增強。
CKC方法避免直接進行有關混合矩陣H的運算，而是試圖從多通道的檢測信號中進行計算求出運動單元的發放序列，而這種計算的主要目的是求出式(7)中一個理想的CX sj，再求出一個結果好的運動單元發放序列，求出CX sj的應能使同一運動單元發放序列的幅值盡可能地提高，這樣就容易觀測出其發放時刻，這就需要保證CX sj應與多通道的信號中同一運動單元發放時刻的信號具有較強的相關性。通過我們自己的實驗驗證發現原有方法提出的求CX sj的方法較繁瑣，而且也不能保證求出的結果是比較滿意的，因而本文提出一種新的方法，通過一個可能同時含有多個運動單元發放的初始序列，這個初始序列的一些最大值所對應的時刻(也即發放脈沖所對應的時刻)往往為一個或多個運動單元的發放序列，利用信號中同一運動單位發放時刻的相似性而采用自組織映射神經網絡具有的無監督聚類的功能來對這些發放時刻的序列進行分類，然后再選出一類進行CX sj的計算。選出的類中所包含的時刻應該全部或是部分屬于同一個運動單位的發放，這樣就能保證所計算出的CX sj能夠較好地計算出一個源發放序列。本文所提出的方法與已有的方法都是努力找出同一類運動單元發放時刻，進而再求出CX sj，與原方法比較，本文方法的計算步驟簡單有效。具體的算法步驟如下：
(1)由C
(2)用式(4)計算γ(n)，并設置，即n0等于γ(n)的中值所對應的時刻，用式(5)計算sn0，找出sn0中的最大值及其所對應的時刻n1，再用式(5)計算sn1，找出sn1中一些(具體個數可依據信號實際情況進行調整)較大值及其所對應的時刻，記這些時刻為φnc={nc1, nc2, …, ncn}。
(3)設置好神經網絡的分類數等具體參數，創建并訓練自組織映射神經網絡，將X(nci), nci∈φnc作為神經網絡的輸入向量，對X(nci)進行分類。
(4)找出含有元素個數最多的一類，記為
(5)設φnv中元素的個數為d, 在中找出(1+m)·d(通常0.5≤m≤5)個與最大值對應的時刻，記為φnt={nt1, nt2, …, ntn}，重新計算，最終得到發放序列為=C
(6)設置γ(d)=0, ?d∈φnv。
(7)重復步驟(2)~(5)Nt次，確保能夠盡量多地重建發放源序列，最后對結果進行歸類處理。
1.3 仿真信號
信號中設置10個運動單元(信號源)[9-10]，隨機產生運動單元的發放序列，其脈沖的平均間隔設置為100個采樣點，Si(k), k=1, 2, …, 200均勻分布在[-10, 10]之間，每一個運動單元發放序列的長度設置為20 000個采樣點，采樣頻率為2 000 Hz。隨機生成零均值混合矩陣H，卷積核的長度P設置為10個采樣點，檢測信號的通道數目M為25，為了增強計算效果，在每一個原檢測通道基礎上再增加9個通道的延遲信號[9-10]，這樣加上擴展的運動單元脈沖序列，運動單元的脈沖序列共增加到190個，檢測通道的數目為250個。每個生成的信號根據需要加上不同信噪比(sigal to noise ratio, SNR)的零均值的高斯白噪聲, 如圖 1所示。
圖1
仿真信號
上圖為用隨機混合矩陣
top panel: one channel of simulated signal generated by the random mixing matrix
2 實驗
2.1 仿真實驗設計
仿真主要研究以下內容：
(1)從初始序列sn1中選取較大值個數不變，研究在較高信噪比和較低信噪比情況下自組織映射神經網絡不同的分類個數對重建出運動單元個數以及重建出的發放脈沖準確度的影響。
(2)在自組織映射神經網絡分類個數不變的情況下，從初始序列sn1中選取不同的較大值個數對重建出運動單元個數以及重建出的發放脈沖準確度的影響。
(3)在自組織映射神經網絡的分類個數以及從初始序列sn1中選取不同的較大值個數不變的情況下，該方法與原有CKC方法在不同的信噪比情況下能夠重建出運動單元的個數和重建出的發放脈沖準確度。
在所有的實驗中，自組織映射神經網絡的訓練次數設置為200次，選擇默認的拓撲結構。
2.2 實驗結果
圖 2展示的是在信噪比為0 dB時，重建出的一個運動單元的發放序列和原始運動單元脈沖發放序列的比較情況，其對幅值進行了歸一化操作，淺色的虛線部分是其中一個信號源的原發放序列，深色部分的為重建出信號源的發放序列。圖 2中可以看出，重建出的序列真實地反映了原始發放序列。
圖2
原發放序列與重建出的發放序列的比較
Figure2.
Comparison between the original pulse train and the reconstructed pulse train
圖 3(a)和(b)展示的是從初始序列sn1選取較大值個數不變以及在信噪比分別為10 dB和-5 dB的情況下，自組織映射神經網絡的分類個數與重建出的運動單元個數及發放脈沖的平均準確度之間的關系。其中選取的初始序列較大值個數都為200，分類個數分別選為3、10、20、30、40、50。由圖可以看出在信噪比為-5 dB時，分類個數為3時，只能重建出3個運動單元序列，并且發放脈沖的平均準確率也偏低；當分類個數在10以上時，重建出的運動單元序列數相差不大，要么為8個，要么為9個，但是難以完全重建出所有10個運動單元，重建出的序列發放脈沖的平均準確度也相差不大。在信噪比為10 dB時，即使分類個數為3的情況下也能重建出所有的運動單元，但是重建出發放的脈沖平均準確率要低于分類個數多的平均準確率，而在分類個數為10以上時，能以相差不大的較高的發放脈沖準確率重建出所有的運動單元。
圖3
重建出的信號源的個數和平均脈沖準確率
(a)不同的分類個數所對應的重建出的信號源個數; (b)不同的分類個數所對應的重建出的發放脈沖平均準確率; (c)選取不同的原始序列較大值個數所對應的重建出的信號源的個數; (d)選取不同的原始序列較大值個數所對應的重建出的發放脈沖平均準確率
Figure3. Number of reconstructed sources and average pulse accuracy(a) number of reconstructed sources versus number of groups clustered by som; (b) average pulse accuracy versus number of groups clustered by som; (c) number of reconstructed sources versus number of highest peaks selected from initial pulse train; (d) average pulse accuracy versus number of highest peaks selected from initial pulse train
圖 3(c)和(d)展示的是在信噪比為0 dB時，經過自組織映射神經網絡分類的個數為6的情況下，從初始序列中選取不同的較大值個數與重建出的運動單元個數及重建發放脈沖的平均準確率之間的關系。較大值的個數分別選為20、50、100、200、300、400、500、600和700。由圖中可以看出，較大值個數選為20的時候，只能重建出6個運動單元，并且重建出的發放的脈沖準確率也較低，較大值個數選為50、100、200、300、400和500時，能重建出所有的運動單元，較大值個數選為600和700時，能重建出9個運動單元，而隨著選取的較大值個數從少到多的變化，重建出的序列發放脈沖的平均準確率也有一個先增加、后減少的過程。
圖 4展示的是不同程度的噪聲對本文所提出的自組織映射的卷積核補償(self-organizing map CKC，SOM CKC)方法和原有卷積核補償(classic CKC)方法各自所能重建出運動單元的個數以及重建序列脈沖發放的平均準確率的影響。信噪比分別設置為-10、-5、0、5和10 dB, 從初始序列中選取的較大值的個數為200，經過自組織映射神經網絡所分類的個數設置為6個。圖 4中可以看出即使信噪比設置為-5 dB時，SOM CKC方法依然能夠能重建出所有的信號源，而重建序列脈沖發放的平均準確率之間相差不是很大，但總體趨勢是信噪比較高時，準確率較高；在信噪比為-10 dB時，該方法能夠重建出4個信號源，但脈沖準確率下降到80%左右。而classic CKC方法在信噪比為-5 dB時能夠重建出8個信號源，在信噪比為-10 dB時難以重建出任何一個信號源。
圖4
兩種不同方法的比較結果
(a)兩種方法在不同信噪比情況下所能重建出的運動單元個數; (b)兩種方法在不同信噪比情況下所能重建出的平均脈沖準確率
Figure4. Comparison results(a) number of sources reconstructed by SOM CKC and classic CKC; (b) average pulse accuracy achieved by SOM CKC and classic CKC
3 討論
本文提出自組織映射神經網絡和已有CKC方法相結合進行多通道信號的信號源的重建和分解工作。該方法步驟簡單有效，經過仿真實驗表明所提出的方法能夠在信號的信噪比即使在-5 dB時，仍能以較高的準確率重建出所有的運動單元。但是當信噪比降到-10 dB時，效果下降，難以重建出所有的運動單元(見圖 4)。
該方法要求經過自組織映射神經網絡分類以后，每一類中包含的發放時刻大部分由同一個運動單元所發放，這樣才能較好地重建出這個運動單元，如果自組織映射神經網絡分類個數設置的較少，分類以后會使同一類中含有大量的非同一個運動單元產生的發放時刻，這樣就不能很好地求出CX sj，勢必會影響到重建的結果，信號的信噪比較高時，能更好地進行分類，而信噪比較低的時候，分類效果會變差，所以無論是重建出的運動單元的個數還是重建出的序列的脈沖準確率來說，高信噪比的結果都會好于低信噪比的結果(見圖 3、4)。如果從初始序列中sn1選取的較大值個數較少，再經過分類以后，每一類中含有的元素個數就會很少，再利用這些元素的平均值求取的CX sj不具和同一脈沖序列的發放時刻的信號具有普遍的相關性，只能和所利用求CX sj的元素
4 結論
本文提出將自組織映射神經網絡和已有的CKC方法進行結合來簡化計算過程并提高信號分解的效率，主要是利用多通道的sEMG信號中，同一運動單元的發放時刻的信號有著一定程度上的相似性，然后利用自組織映射神經網絡的無監督學習的聚類功能，對初始序列中的一些較大值對應的時刻進行分類，進而可求出CX sj，這樣最終就可以重建信號源的發放序列。由隨機混合矩陣產生的仿真信號實驗表明，該方法能在信號信噪比較低的情況下，以較高的準確率重建出幾乎所有信號源的發放序列，證明了本文所提出的方法的有效性。
引言
生物醫學信號往往包含著復雜的信息，通常由多個信號源的信號綜合所組成，比如神經元和神經纖維。這些信號源通常被看作是統計獨立的，或者至少是弱相關的，但是它們的混合過程實際上是不知道的，因此需要對所獲得的這類信號進行盲分解。
過去一段時間以來針對肌電(electromyogram, EMG)信號的研究是很活躍的，因為無論其在肌肉康復治療中還是在臨床診斷中都有著重要的現實意義，研究者們還通過測量EMG信號研究神經肌肉的控制機理以及在其它方面進行應用，例如假肢控制[1-2]，手語識別[3]等。EMG信號中的動作電位(action potential, AP)通常由多個運動單位(motor unit, MU)所產生的，其被看作是隨機干擾的非平穩信號，具有復雜的特性，尤其是表面肌電(surface electromyogram, sEMG)信號的復雜程度更要高一些，不但其測量信號的電極與所接觸皮膚的接觸特性不良好，而且還有各種噪聲的干擾，因此加劇了對sEMG信號分解和提取臨床相關信息的困難。
近年來，高密度陣列式的sEMG電極被用來進行sEMG信號的檢測，針對這種檢測方式，一些sEMG信號處理的方法被陸續提出，包括模式識別的方法和盲源分離的方法。Kleine等[4]指出二維空間濾波器在sEMG信號分解中對檢測發放精度的重要性。Gazzoni等[5]采用分類技術和模板匹配方法對信號進行分解研究。Wood等[6]通過有限元分析方法來研究運動單元的深度和皮膚表面電位分布之間的關系。Garcia等[7]把sEMG信號的產生看作是運動單元動作電位(motor unit action potential, MUAP)序列的瞬時混合過程。Gligorijevic等[8]先對運動單元電位波形進行聚類，然后再用一種分解技術對信號進行分解。Holobar等[9-13]提出了卷積核補償(gradient convolution kernel compensation, CKC)方法和梯度優化的CKC方法來處理卷積混合的運動單元發放源信號，這種方法不需要計算未知的混合矩陣并且對發放源信號直接進行計算。CKC方法被證明是很有效的，如果sEMG信號的質量較好，可以同時重建出30多個發放源[10]。
本文在現有CKC方法的理論研究和實踐應用的基礎之上，嘗試提出一種新的方法，該方法將結合自組織映射神經網絡對信號進行處理，以期對現有的CKC方法的計算過程進行進一步的簡化并獲得更好的處理結果。
1 原理與方法
1.1 數據模型和卷積核補償
對于一個線性時不變(linear time invariable, LTI)的多輸入多輸出(multi-input multi-output, MIMO)系統，CKC的數學模型可以用以下方程表示[9-10]：
| $ \boldsymbol{X}\text{ }(n)=\boldsymbol{H}\bar{s}\left(n \right)+\boldsymbol{e}\left(n \right), $ |
其中X(n)=[x1(n)，…, xM(n)]T為卷積后的多通道檢測信號向量，xj(n)為第j個通道中的第n個采樣值，e(n)=[e1(n), …, eM(n)]T為零均值白噪聲向量，
| $ \boldsymbol{H}=\left(\begin{align} & {{h}_{11}}\left(0 \right)\cdots {{h}_{11}}\left(P-1 \right)\ \ \ {{h}_{12}}\left(0 \right)\cdots {{h}_{12}}\left(P-1 \right)\cdots {{h}_{1N}}\left(0 \right)\cdots {{h}_{1N}}\left(P-1 \right)\\ & {{h}_{21}}\left(0 \right)\cdots {{h}_{21}}\left(P-1 \right)\ \ \ {{h}_{22}}\left(0 \right)\cdots {{h}_{22}}\left(P-1 \right)\cdots {{h}_{2N}}\left(0 \right)\cdots {{h}_{2N}}\left(P-1 \right)\\ & \ \ \ \vdots \ \ \ \ \ \cdots \ \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \ \ \cdots \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \cdots \ \ \ \ \vdots \ \ \ \ \ \ \ \cdots \ \ \ \ \ \ \ \ \ \ \vdots \\ & {{h}_{M1}}\left(0 \right)\cdots {{h}_{M1}}\left(P-1 \right)\ \ \ {{h}_{M2}}\left(0 \right)\cdots {{h}_{M2}}\left(P-1 \right)\cdots {{h}_{MN}}\left(0 \right)\cdots {{h}_{MN}}\left(P-1 \right)\\ \end{align} \right), $ |
其中hij是在第i個通道中第j個信號源的源脈沖序列的響應，si(n)為第i個信號源的發放序列，可以用以下脈沖序列的模型表示：
| $ {{s}_{i}}\left(n \right)=\sum\limits_{-\infty }^{\infty }{\delta \left[n-{{S}_{i}}\left(k \right)\right]}, $ |
其中δ(·)為狄拉克沖激函數，Si(k)為第i個信號源的第k個發放時刻。
為了改善通道數和未知信號源的數目之間的比例，以提高分解效果，在模型(1)中，對每一個X(n)分量進行k-1個延遲擴展，即
| $ \begin{align} & \gamma \left(n \right)={{{\bar{\boldsymbol{X}}}}^{\boldsymbol{T}}}\left(n \right)\boldsymbol{C}_{\bar{\boldsymbol{X}}\bar{\boldsymbol{X}}}^{-\ 1}\bar{\boldsymbol{X}}\left(n \right)=\bar{s}\left(n \right){{\boldsymbol{H}}^{T}}{{\left({{\boldsymbol{H}}^{T}} \right)}^{-1}}\boldsymbol{C}_{\bar{s}\bar{s}}^{-1}{{\boldsymbol{H}}^{-1}} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \boldsymbol{H}\ \bar{s}\left(n \right)={{{\bar{s}}}^{T}}\left(n \right)\boldsymbol{C}_{\bar{s}\bar{s}}^{-1}\bar{s}\left(n \right), \\ \end{align} $ |
其中C
如果式(4)中右側第一項用
| $ {{s}_{{{n}_{0}}}}\left(n \right)={{{\bar{\boldsymbol{X}}}}^{T}}\left({{n}_{0}} \right)\boldsymbol{C}_{\bar{\boldsymbol{X}}\bar{\boldsymbol{X}}}^{-\ 1}\bar{\boldsymbol{X}}\left(n \right), $ |
已經證明如果第j個信號源在時刻n0發放，sn0(n)可表示為對第j個信號源發放序列的估計[9]，
| $ {{s}_{{{n}_{0}}}}\left(n \right)\approx {{{\bar{s}}}_{j}}\left(n \right), $ |
這樣可以通過式(5)計算出一個信號源的脈沖發放序列，以進行信號的分解。為了改善sn0(n)的計算質量，式(5)中的
| $ {{\boldsymbol{C}}_{{{\boldsymbol{X}}_{sj}}}}=\frac{1}{card\left({{\varphi }_{n}} \right)}\sum{\bar{\boldsymbol{X}}\left({{\varphi }_{n}} \right)}, $ |
其中φn={n0, n1, …, nm}，card(φn)表示φn中所包含元素的數目。φn中所有的時刻應該被同一個信號源的脈沖序列的發放時刻所包括，或是絕大部分被包括，這樣計算出的效果會更好。
1.2 基于自組織映射神經網絡的卷積核補償方法
自組織映射神經網絡[14-16]是基于無監督學習方法的神經網絡的一種重要類型。自組織映射網絡理論最早是由芬蘭赫爾辛基理工大學的Kohohen于1981年提出的。此后，伴隨著神經網絡在20世紀80年代中后期的迅速發展，自組織映射理論及其應用也有了長足的進步。
自組織影射神經網絡的訓練過程是無監督的聚類過程, 其通過對輸入樣本的反復學習, 使連接權向量的空間分布與輸入樣本的概率分布趨于一致, 進而將多維的任意輸入模式在輸出層映射成一維、二維或多維離散圖形, 并且圖形保持輸入的拓撲結構不變。自組織映射網絡的學習算法抽象意義上包括：競爭、合作和更新三個階段。
(1)競爭：對于每一個輸入樣本, 網絡中的神經元對它們各自的判別函數值進行計算, 這個判別函數值對神經元之間的競爭提供了基礎, 最終具有函數最大值的神經元將成為競爭的勝利者。
(2)合作：獲勝的神經元決定了興奮神經元的拓撲鄰域空間位置, 進而為鄰域神經元提供合作的基礎。
(3)更新：興奮神經元通過對它們各自權值的適當調節以增強該輸入模式的判別函數值, 使以后相似模式輸入的響應增強。
CKC方法避免直接進行有關混合矩陣H的運算，而是試圖從多通道的檢測信號中進行計算求出運動單元的發放序列，而這種計算的主要目的是求出式(7)中一個理想的CX sj，再求出一個結果好的運動單元發放序列，求出CX sj的應能使同一運動單元發放序列的幅值盡可能地提高，這樣就容易觀測出其發放時刻，這就需要保證CX sj應與多通道的信號中同一運動單元發放時刻的信號具有較強的相關性。通過我們自己的實驗驗證發現原有方法提出的求CX sj的方法較繁瑣，而且也不能保證求出的結果是比較滿意的，因而本文提出一種新的方法，通過一個可能同時含有多個運動單元發放的初始序列，這個初始序列的一些最大值所對應的時刻(也即發放脈沖所對應的時刻)往往為一個或多個運動單元的發放序列，利用信號中同一運動單位發放時刻的相似性而采用自組織映射神經網絡具有的無監督聚類的功能來對這些發放時刻的序列進行分類，然后再選出一類進行CX sj的計算。選出的類中所包含的時刻應該全部或是部分屬于同一個運動單位的發放，這樣就能保證所計算出的CX sj能夠較好地計算出一個源發放序列。本文所提出的方法與已有的方法都是努力找出同一類運動單元發放時刻，進而再求出CX sj，與原方法比較，本文方法的計算步驟簡單有效。具體的算法步驟如下：
(1)由C
(2)用式(4)計算γ(n)，并設置，即n0等于γ(n)的中值所對應的時刻，用式(5)計算sn0，找出sn0中的最大值及其所對應的時刻n1，再用式(5)計算sn1，找出sn1中一些(具體個數可依據信號實際情況進行調整)較大值及其所對應的時刻，記這些時刻為φnc={nc1, nc2, …, ncn}。
(3)設置好神經網絡的分類數等具體參數，創建并訓練自組織映射神經網絡，將X(nci), nci∈φnc作為神經網絡的輸入向量，對X(nci)進行分類。
(4)找出含有元素個數最多的一類，記為
(5)設φnv中元素的個數為d, 在中找出(1+m)·d(通常0.5≤m≤5)個與最大值對應的時刻，記為φnt={nt1, nt2, …, ntn}，重新計算，最終得到發放序列為=C
(6)設置γ(d)=0, ?d∈φnv。
(7)重復步驟(2)~(5)Nt次，確保能夠盡量多地重建發放源序列，最后對結果進行歸類處理。
1.3 仿真信號
信號中設置10個運動單元(信號源)[9-10]，隨機產生運動單元的發放序列，其脈沖的平均間隔設置為100個采樣點，Si(k), k=1, 2, …, 200均勻分布在[-10, 10]之間，每一個運動單元發放序列的長度設置為20 000個采樣點，采樣頻率為2 000 Hz。隨機生成零均值混合矩陣H，卷積核的長度P設置為10個采樣點，檢測信號的通道數目M為25，為了增強計算效果，在每一個原檢測通道基礎上再增加9個通道的延遲信號[9-10]，這樣加上擴展的運動單元脈沖序列，運動單元的脈沖序列共增加到190個，檢測通道的數目為250個。每個生成的信號根據需要加上不同信噪比(sigal to noise ratio, SNR)的零均值的高斯白噪聲, 如圖 1所示。
圖1
仿真信號
上圖為用隨機混合矩陣
top panel: one channel of simulated signal generated by the random mixing matrix
2 實驗
2.1 仿真實驗設計
仿真主要研究以下內容：
(1)從初始序列sn1中選取較大值個數不變，研究在較高信噪比和較低信噪比情況下自組織映射神經網絡不同的分類個數對重建出運動單元個數以及重建出的發放脈沖準確度的影響。
(2)在自組織映射神經網絡分類個數不變的情況下，從初始序列sn1中選取不同的較大值個數對重建出運動單元個數以及重建出的發放脈沖準確度的影響。
(3)在自組織映射神經網絡的分類個數以及從初始序列sn1中選取不同的較大值個數不變的情況下，該方法與原有CKC方法在不同的信噪比情況下能夠重建出運動單元的個數和重建出的發放脈沖準確度。
在所有的實驗中，自組織映射神經網絡的訓練次數設置為200次，選擇默認的拓撲結構。
2.2 實驗結果
圖 2展示的是在信噪比為0 dB時，重建出的一個運動單元的發放序列和原始運動單元脈沖發放序列的比較情況，其對幅值進行了歸一化操作，淺色的虛線部分是其中一個信號源的原發放序列，深色部分的為重建出信號源的發放序列。圖 2中可以看出，重建出的序列真實地反映了原始發放序列。
圖2
原發放序列與重建出的發放序列的比較
Figure2.
Comparison between the original pulse train and the reconstructed pulse train
圖 3(a)和(b)展示的是從初始序列sn1選取較大值個數不變以及在信噪比分別為10 dB和-5 dB的情況下，自組織映射神經網絡的分類個數與重建出的運動單元個數及發放脈沖的平均準確度之間的關系。其中選取的初始序列較大值個數都為200，分類個數分別選為3、10、20、30、40、50。由圖可以看出在信噪比為-5 dB時，分類個數為3時，只能重建出3個運動單元序列，并且發放脈沖的平均準確率也偏低；當分類個數在10以上時，重建出的運動單元序列數相差不大，要么為8個，要么為9個，但是難以完全重建出所有10個運動單元，重建出的序列發放脈沖的平均準確度也相差不大。在信噪比為10 dB時，即使分類個數為3的情況下也能重建出所有的運動單元，但是重建出發放的脈沖平均準確率要低于分類個數多的平均準確率，而在分類個數為10以上時，能以相差不大的較高的發放脈沖準確率重建出所有的運動單元。
圖3
重建出的信號源的個數和平均脈沖準確率
(a)不同的分類個數所對應的重建出的信號源個數; (b)不同的分類個數所對應的重建出的發放脈沖平均準確率; (c)選取不同的原始序列較大值個數所對應的重建出的信號源的個數; (d)選取不同的原始序列較大值個數所對應的重建出的發放脈沖平均準確率
Figure3. Number of reconstructed sources and average pulse accuracy(a) number of reconstructed sources versus number of groups clustered by som; (b) average pulse accuracy versus number of groups clustered by som; (c) number of reconstructed sources versus number of highest peaks selected from initial pulse train; (d) average pulse accuracy versus number of highest peaks selected from initial pulse train
圖 3(c)和(d)展示的是在信噪比為0 dB時，經過自組織映射神經網絡分類的個數為6的情況下，從初始序列中選取不同的較大值個數與重建出的運動單元個數及重建發放脈沖的平均準確率之間的關系。較大值的個數分別選為20、50、100、200、300、400、500、600和700。由圖中可以看出，較大值個數選為20的時候，只能重建出6個運動單元，并且重建出的發放的脈沖準確率也較低，較大值個數選為50、100、200、300、400和500時，能重建出所有的運動單元，較大值個數選為600和700時，能重建出9個運動單元，而隨著選取的較大值個數從少到多的變化，重建出的序列發放脈沖的平均準確率也有一個先增加、后減少的過程。
圖 4展示的是不同程度的噪聲對本文所提出的自組織映射的卷積核補償(self-organizing map CKC，SOM CKC)方法和原有卷積核補償(classic CKC)方法各自所能重建出運動單元的個數以及重建序列脈沖發放的平均準確率的影響。信噪比分別設置為-10、-5、0、5和10 dB, 從初始序列中選取的較大值的個數為200，經過自組織映射神經網絡所分類的個數設置為6個。圖 4中可以看出即使信噪比設置為-5 dB時，SOM CKC方法依然能夠能重建出所有的信號源，而重建序列脈沖發放的平均準確率之間相差不是很大，但總體趨勢是信噪比較高時，準確率較高；在信噪比為-10 dB時，該方法能夠重建出4個信號源，但脈沖準確率下降到80%左右。而classic CKC方法在信噪比為-5 dB時能夠重建出8個信號源，在信噪比為-10 dB時難以重建出任何一個信號源。
圖4
兩種不同方法的比較結果
(a)兩種方法在不同信噪比情況下所能重建出的運動單元個數; (b)兩種方法在不同信噪比情況下所能重建出的平均脈沖準確率
Figure4. Comparison results(a) number of sources reconstructed by SOM CKC and classic CKC; (b) average pulse accuracy achieved by SOM CKC and classic CKC
3 討論
本文提出自組織映射神經網絡和已有CKC方法相結合進行多通道信號的信號源的重建和分解工作。該方法步驟簡單有效，經過仿真實驗表明所提出的方法能夠在信號的信噪比即使在-5 dB時，仍能以較高的準確率重建出所有的運動單元。但是當信噪比降到-10 dB時，效果下降，難以重建出所有的運動單元(見圖 4)。
該方法要求經過自組織映射神經網絡分類以后，每一類中包含的發放時刻大部分由同一個運動單元所發放，這樣才能較好地重建出這個運動單元，如果自組織映射神經網絡分類個數設置的較少，分類以后會使同一類中含有大量的非同一個運動單元產生的發放時刻，這樣就不能很好地求出CX sj，勢必會影響到重建的結果，信號的信噪比較高時，能更好地進行分類，而信噪比較低的時候，分類效果會變差，所以無論是重建出的運動單元的個數還是重建出的序列的脈沖準確率來說，高信噪比的結果都會好于低信噪比的結果(見圖 3、4)。如果從初始序列中sn1選取的較大值個數較少，再經過分類以后，每一類中含有的元素個數就會很少，再利用這些元素的平均值求取的CX sj不具和同一脈沖序列的發放時刻的信號具有普遍的相關性，只能和所利用求CX sj的元素
4 結論
本文提出將自組織映射神經網絡和已有的CKC方法進行結合來簡化計算過程并提高信號分解的效率，主要是利用多通道的sEMG信號中，同一運動單元的發放時刻的信號有著一定程度上的相似性，然后利用自組織映射神經網絡的無監督學習的聚類功能，對初始序列中的一些較大值對應的時刻進行分類，進而可求出CX sj，這樣最終就可以重建信號源的發放序列。由隨機混合矩陣產生的仿真信號實驗表明，該方法能在信號信噪比較低的情況下，以較高的準確率重建出幾乎所有信號源的發放序列，證明了本文所提出的方法的有效性。
