用于等離子體相干模式自動識別的譜聚類算法實現(xiàn)

趙子博 莊革 謝錦林 渠承明 強子薇

(中國科學(xué)技術(shù)大學(xué)核科學(xué)技術(shù)學(xué)院,合肥 230026)

高約束模式對改善等離子體約束有著重要意義,但目前主要依賴人工進行模式識別,其效率低、成本高,導(dǎo)致核聚變裝置中大量的診斷數(shù)據(jù)沒有得到充分分析.為了解決這個問題,本文將機器學(xué)習(xí)中的譜聚類算法應(yīng)用到EAST 托卡馬克裝置上的電子回旋輻射成像、一維診斷系統(tǒng)電子回旋輻射計、磁探針、軟X 射線和快輻射等不同診斷系統(tǒng)的數(shù)據(jù)上,在時域及頻域上識別出了鋸齒模,驗證了譜聚類方法的遷移性及準確性,解決了監(jiān)督學(xué)習(xí)在數(shù)據(jù)處理上遷移性差以及需要依賴大量標簽數(shù)據(jù)的問題.此外,本文實現(xiàn)了特定模式的篩選;最后利用電子回旋輻射成像及磁探針數(shù)據(jù)發(fā)現(xiàn)了一種可能的新模式,為新模式探索提供了一種新思路.

1 引言

自1982 年ASDEX 裝置第一次獲得了高約束模式[1]以來,托卡馬克裝置的能量約束時間不斷提高.由于高約束模式對等離子體具有良好的約束性能,因此被認為是最有可能實現(xiàn)核聚變的運行模式.但同時在高約束情況下存在各種不穩(wěn)定模式,比如伴隨著臺基區(qū)等離子體約束性能的周期性下降的邊界局域模[2];伴隨著等離子體芯部密度和溫度的周期性耗散的鋸齒模[3].為了優(yōu)化托卡馬克的設(shè)計以改善等離子體約束,必須要對等離子體中的模式進行識別.

核聚變裝置已經(jīng)積累了大量的診斷數(shù)據(jù),例如EAST 托卡馬克裝置上的電子回旋輻射成像(ECEI)系統(tǒng)自2012 年以來,已經(jīng)采集超過7000炮,每一炮的數(shù)據(jù)大小約為7.6 GB,總數(shù)據(jù)量已超過40 TB[1].此外,診斷系統(tǒng)繁多,比如一維診斷系統(tǒng)電子回旋輻射計(ECE)[4]和 ECEI[5?9]等,因此總數(shù)據(jù)量巨大.然而目前主要依賴人工進行模式區(qū)分,該方法效率低、成本高,無法滿足實際需求,導(dǎo)致大量的診斷數(shù)據(jù)沒有得到充分分析.因此尋找高效模式識別的方法十分重要.

近年來,機器學(xué)習(xí)在核聚變領(lǐng)域已經(jīng)有了廣泛而深入的應(yīng)用.機器學(xué)習(xí)分為監(jiān)督學(xué)習(xí)與無監(jiān)督學(xué)習(xí),目前在可控核聚變領(lǐng)域應(yīng)用最多的是監(jiān)督學(xué)習(xí)[10?12].但是監(jiān)督學(xué)習(xí)在處理數(shù)據(jù)方面有很大的缺陷,一方面監(jiān)督學(xué)習(xí)需要大量帶有標簽的數(shù)據(jù)來進行訓(xùn)練,而目前大多數(shù)診斷數(shù)據(jù)的標簽尚未完善;另一方面監(jiān)督學(xué)習(xí)的遷移性很差,對于不同裝置或者同一裝置上的不同診斷系統(tǒng),甚至同一裝置上的同一個診斷系統(tǒng)不同條件下的數(shù)據(jù)都需要重新訓(xùn)練模型,效率低、適用性差.而無監(jiān)督學(xué)習(xí)不存在此缺點,因此本文采用無監(jiān)督學(xué)習(xí)的方法.

目前無監(jiān)督學(xué)習(xí)以聚類算法為主,傳統(tǒng)的聚類算法有K均值算法[13]、層次聚類算法以及密度聚類算法.其中,K均值算法對初值敏感,僅適用于凸形簇[14];密度聚類時間復(fù)雜度高、效率低,參數(shù)選取缺乏理論性;層次聚類算法不能更正錯誤的決策且偏好凸形簇,且對簇的大小有一定的要求.相比之下,基于圖論的譜聚類算法首先具備遷移性強的特點;其次,該算法對樣本分布的適應(yīng)性強[1],可以用來識別各種形狀的簇;最后,譜聚類算法對初值不敏感[15,16],且應(yīng)用該算法可以高效且準確地處理數(shù)據(jù)[1].

綜上,本文主要以譜聚類算法為基本算法,對EAST 托卡馬克裝置上不同診斷系統(tǒng)包括ECEI、ECE、磁探針、軟X 射線(SXR)和快輻射,密度診斷上的數(shù)據(jù)進行自動處理,實現(xiàn)自動尋找模式的目的,為新模式探索提供一種新思路.同時還能實現(xiàn)特定模式的篩選,大幅減少研究人員用于數(shù)據(jù)處理上的時間.

2 譜聚類算法實現(xiàn)

2.1 類別數(shù)已知情況下算法原理

聚類算法的任務(wù)是基于數(shù)據(jù)間的關(guān)系將不同的樣本劃分成多個不相交子集.譜聚類算法將樣本看成空間的點,每兩個點之間用一條被賦予權(quán)值的邊連接.每條邊的權(quán)值表示樣本之間的相似度,權(quán)值越大,相似度越強.通過對該圖的劃分,使各個子圖內(nèi)部邊權(quán)的和越大越好,不同子圖間邊權(quán)的和越小越好,進而實現(xiàn)聚類的目的[1].假設(shè)N個樣本,K個類別,定義損失函數(shù)為[9]

其中wij表示第i個點與第j個點之間的相似度,Aj表示第j個類別,vi表示第i個樣本,di表示所有樣本與第i個樣本之間的相似度的總和,vol(Al) 表示屬于該類別的樣本與所有樣本之間相似度的總和.

損失函數(shù)最小對應(yīng)于最合理的分類標準.但是由于尋找損失函數(shù)最小值是一個無法在多項式時間內(nèi)計算求解的問題(NP 難問題),因此要尋找一種近似算法使其在有限的計算資源和時間下可以求解.首先,將損失函數(shù)改寫為矩陣形式:

進一步,定義矩陣:

定義

求解損失函數(shù)最小值,即為求解 tr(PTV P)在PTP為單位矩陣條件下的最大值.對于這類條件極值問題,類似主成分分析,可以運用拉格朗日乘子法進行求解,即

其中β是待定常數(shù).通過(11)式和(12)可以得到

求解該特征方程可以得到N個特征值,即β的N個取值.其中最大的K個特征值的和就是損失函數(shù)的最大值,即

將選取的特征值所對應(yīng)的特征向量組合構(gòu)成P矩陣.令P矩陣的K個列向量組成K維子空間.根據(jù)P矩陣的定義,P矩陣包含著分類信息,每個行向量與一個樣本對應(yīng),同類樣本趨于K維子空間的一個軸分布,不同類樣本會在不同軸上分布.計算每個樣本(P矩陣的每個行向量)到每個類別聚類中心的距離,用 (vj,ci)dist來表示第j個樣本與第i個聚類中心間的距離,

其中,vj表示第j個樣本,ci表示第i個類別的聚類中心,n表示屬于第i類的樣本數(shù).將該樣本劃分到與它距離最小的聚類中心所代表的類別中,所有樣本分配完畢后,重新計算聚類中心.如果聚類中心發(fā)生變化,重新分配樣本直到聚類中心不再發(fā)生變化為止.值得注意的是,以上過程體現(xiàn)了譜聚類算法的兩個優(yōu)勢,其中一個優(yōu)勢是將原本復(fù)雜的樣本結(jié)構(gòu)轉(zhuǎn)換成了簡單的分布(同類樣本在坐標系中的一條直線上分布,不同類樣本在不同直線上分布),便于分類,因此保證了算法的準確性;另一個優(yōu)勢是將原本高維的樣本數(shù)據(jù)進行了降維(P矩陣的每個行向量的維度是類別數(shù)),提高了算法的效率.

2.2 類別數(shù)未知情況下算法原理

在實際工作中,多數(shù)情況下類別數(shù)不能事先確定,因此需要一個方法自動確定類別數(shù).定義第j個樣本與第i個聚類中心間的距離為

其中

ε為一個很小的數(shù),取為eps,即ε=2.2204×10–16;γ用來控制簇的粗細,γ越小,分類標準越高,通常取0.01.

按上述定義距離的方法,可以使得同種類別(徑向分布)樣本之間的距離小于樣本到坐標原點的距離,不同種類樣本之間的距離大于樣本到坐標原點的距離[1].因此,在最開始進行分類的時候,可以假定樣本會被分成3 類,其中2 個類別的聚類中心分別由樣本間相似度最低的2 個樣本定義(確保2 個樣本不是同一類);第3 個類別的聚類中心為坐標原點.之后,計算每個樣本分別到3 個聚類中心的距離,將樣本劃分到與它距離最小的聚類中心所代表的類別中.所有樣本劃分完成后,更新聚類中心,循環(huán)迭代,直到聚類中心不再變化為止.完成后,如果原點所代表的類別中沒有樣本,則分類完成,類別數(shù)為兩類;如果有樣本,說明類別數(shù)不止兩類,需要將類別數(shù)調(diào)整為4,重復(fù)上述過程,直到K+1 類時,原點所代表的類別里無樣本,則分類完成,類別數(shù)為K.

2.3 算法實現(xiàn)流程

根據(jù)2.1 節(jié)和2.2 節(jié)的分析,譜聚類算法的基本流程為圖1 所示.

圖1 譜聚類算法流程圖Fig.1.Flow chart of spectral clustering algorithm.

3 等離子體相干模式識別及特定模式篩選

3.1 多維診斷數(shù)據(jù)空間聚類

為了對等離子體模式的空間特征進行更好的研究,發(fā)展出了大量的多維診斷系統(tǒng),比如ECEI和SXR 成像陣列等.多維診斷系統(tǒng)可以給出空間各點的信息,因此可以利用空間聚類的方法對空間各點進行分類,每一種類別對應(yīng)一種模式,以此來尋找其中的模式.本文以ECEI 診斷數(shù)據(jù)識別為例進行說明.

EAST 托卡馬克裝置上的ECEI 診斷系統(tǒng)有24 行、16 列獨立的數(shù)據(jù)信道[17?19].每個數(shù)據(jù)信道對應(yīng)一個樣本,總共有384 個樣本.每一個樣本是一個時間序列,兩個樣本之間的相似度為

采取這種定義方式是因為同一類別(模式)的表現(xiàn)形式是數(shù)據(jù)隨時間的變化規(guī)律相同,但幅值可以不同.余弦距離關(guān)注方向上的變化,不關(guān)注幅值,采取余弦距離定義相似度正好符合這個表現(xiàn)形式.之后按圖1 所示的操作流程對ECEI 數(shù)據(jù)每隔0.1 s進行一次聚類識別.在42987 炮1.3—9.4 s 內(nèi)識別出的模式如圖2 所示,每個方格對應(yīng)一個ECEI 的信道,總共384 個;白色與黑色各代表1 種類別.

圖2 聚類識別分類結(jié)果Fig.2.Cluster recognition classification results.

為了證明所識別出來的確實是一種模式,現(xiàn)在以4.0—4.1 s 為例對其進行物理上的一些分析.首先從白色區(qū)域選出一個信道A (第12 行,第9 列),再從黑色區(qū)域選出一個信道C (第2 行,第4 列),最后在白色區(qū)域與黑色區(qū)域交界處選出一個信道B (第20 行,第9 列)畫出時序圖,如圖3 所示,其中,δTe/Te(Te?〈Te〉)/〈Te〉,Te代表對應(yīng)時刻的電子溫度,〈Te〉為4—4.1 s 內(nèi)電子溫度的平均值.可知信道A 與信道C 信號明顯分為爬升期、先兆振蕩期和快速崩塌期3 個階段,符合鋸齒模[5,20,21]的基本特征,因此可以判斷出A 通道信號為正鋸齒,C 通道信號為反鋸齒;交界處B 通道信號溫度保持不變,為反轉(zhuǎn)半徑位置.正鋸齒和反鋸齒的同時存在可以視為判斷鋸齒不穩(wěn)定性的簡單依據(jù)[5].

圖3 各信道的信號時序圖 (a) 信道A;(b) 信道B;(c) 信道CFig.3.Signal timing diagram of the different channel:(a) Channel A;(b) channel B;(c) channel C.

圖4 給出了ECEI 觀測到的鋸齒不穩(wěn)定性演化過程圖,標號(1)—(8) 依次對應(yīng)8 個時刻點;圖4(b)黑色、紅色、藍色曲線分別代表圖4(a)各圖對應(yīng)顏色點處的時序圖.從圖4(b)可以明顯看出:最開始鋸齒崩塌結(jié)束,冷磁島占據(jù)整個q=1 面;隨著等離子體加熱,鋸齒爬升,芯部電子溫度緩慢提高,之后重聯(lián)發(fā)生,芯部熱量向外輸運;最后鋸齒崩塌,冷磁島重新占據(jù)整個q=1 面,符合鋸齒模的演化過程.從圖4(a)可以看出,整個演化圖的空間結(jié)構(gòu)與利用譜聚類的方法識別出的模式空間結(jié)構(gòu)基本一致,說明識別出的白色區(qū)域?qū)?yīng)反轉(zhuǎn)半徑以內(nèi)的區(qū)域,為正鋸齒;黑色區(qū)域?qū)?yīng)反轉(zhuǎn)半徑與混合半徑之間的區(qū)域,為反鋸齒,證明了譜聚類方法的可靠性.

圖4 (a)鋸齒?？臻g結(jié)構(gòu)隨時間的演化過程;(b) 黑色、紅色、藍色曲線分別代表圖4(a)各圖對應(yīng)顏色點處的時序圖Fig.4.(a) Evolution of the space structure of sawtooth mode with time;(b) the black,red,and blue curves respectively represent the timing diagrams at the corresponding color points of each panel in Fig.4(a).

用查準率P和查全率R來衡量算法的準確性,定義為

其中,TP表示真正例,FN表示假反例,FP表示假正例.TP指真實情況和識別結(jié)果均為正例;FP指識別結(jié)果為正例,但真實結(jié)果為反例;FN指識別結(jié)果為反例,但真實結(jié)果為正例.對所有的識別結(jié)果進行統(tǒng)計,在42987 炮的38400 個時間片段(384 個信道,每隔0.1 s 聚類一次,數(shù)據(jù)采集時間10 s)內(nèi),聚類的結(jié)果顯示共有13041 個時間片段被識別為正鋸齒,共有9558 個時間片段被識別為反鋸齒.通過測量的信號時序圖,可以判定實際存在正鋸齒的時間片段共有12555 個,實際存在反鋸齒的時間片段共有9234 個.此外,可以判定在聚類算法識別判定為正鋸齒的13041 個時間片段中共有12150 個片段是真實的正鋸齒,在聚類算法識別判定為反鋸齒的9558 個時間片段中共有8829個片段是真實的反鋸齒.根據(jù)查準率與查全率的定義,可以計算得到正鋸齒的查全率為96.8%,查準率為93.2%;反鋸齒的查全率為95.6%,查準率為92.4%.以上結(jié)果表明,譜聚類算法在識別準確性上表現(xiàn)良好.

3.2 一維診斷數(shù)據(jù)時間聚類

在核聚變裝置上除了多維診斷系統(tǒng),還存在大量一維診斷系統(tǒng),包括ECE、磁探針、弦積分密度測量、SXR 以及快輻射等.一維診斷數(shù)據(jù)通常反映空間單點或者單通道的信息,相比多維診斷系統(tǒng),一維診斷系統(tǒng)可供分類的信息更少,可以利用時間聚類來自動識別其中的相干模式.下面以ECE 診斷系統(tǒng)為例進行具體說明.

對ECE 診斷的時序信號進行傅里葉變換,得到各個時間點的頻率信息,每個時間點的頻率序列對應(yīng)一個樣本,同種模式的表現(xiàn)特征是頻率序列強度相似.各個樣本的相似度用樣本之間的指數(shù)距離表示為

其中σ為人為規(guī)定的參數(shù),用來控制樣本間的相似度,本文中σ2=1000.

對數(shù)據(jù)每隔0.1 s 識別一次,在50015 炮的1.8—9.5 s 內(nèi)發(fā)現(xiàn)了一種模式,圖5 為該模式的頻譜圖.通過頻譜圖,可以發(fā)現(xiàn)有展寬非常大的破裂,符合鋸齒模的特征,認定識別出的模式為鋸齒模.對于SXR、快輻射以及ECE 三種診斷數(shù)據(jù)仿照前述操作進行時間聚類,對50015—50115 炮的24000 個時間片段的識別結(jié)果進行統(tǒng)計,其中6700 個時間片段被識別為鋸齒模.根據(jù)頻譜圖,可以判定有6555 個真實存在鋸齒模的時間片段.同時,根據(jù)頻譜圖也可以判定出在被聚類算法識別出的6700 個鋸齒模片段中,有6057 個是真實的鋸齒模片段.因此,根據(jù)(20)式和(21)式對查準率和查全率的定義,可以計算得到P=90.4%,R=92.4%.

圖5 對于50015 炮,模式頻率特征Fig.5.Mode frequency characteristics for shot 50015.

為了實現(xiàn)自動篩選模式,將聚類識別找到的在50015 炮1.8—9.5 s 內(nèi)的鋸齒模的典型信號,即聚類中心提取出來加到時間聚類的樣本里,并作為初始聚類中心.在識別過程中,與該序列分在一類的便是該種模式,以此達到篩選特定模式的功能.對42987—50180 炮內(nèi)的480000 個時間片段進行篩選識別,其中10730 個時間片段被識別為鋸齒模.根據(jù)頻譜圖,可以判定有10719 個真實存在鋸齒模的時間片段;同時,根據(jù)頻譜圖也可以判定出在被聚類算法識別出的10730 個鋸齒模片段中有10708 個是真實的鋸齒模片段.因此,根據(jù)(20)式和(21)式,可以計算得到查準率為99.8%,查全率為99.9%.證明譜聚類算法在模式篩選上的表現(xiàn)非常好,可以實現(xiàn)篩選特定模式的功能,大幅減少研究人員的時間.

3.3 新模式探索

仿照前述鋸齒模的識別方法,利用ECEI 數(shù)據(jù)在64960 炮3.3—3.6 s 內(nèi)發(fā)現(xiàn)了一種模式,見圖6.圖6(a)中每個方格對應(yīng)一個信道,總共有384 個;對384 個樣本進行空間聚類,識別出了一種模式,用淺藍色方格表示.圖6(b)是通過時域圖及頻譜圖判斷出的該模式實際出現(xiàn)的位置,藍色方格區(qū)代表模式出現(xiàn)的地方.兩者對比,發(fā)現(xiàn)聚類識別出的結(jié)果與模式實際出現(xiàn)的位置基本吻合.圖6(c)反映該模式在托卡馬克裝置中的實際位置,其中黑色方格區(qū)即為圖6(b)在托卡馬克裝置中的實際位置.EFIT 代表磁面,信號代表模式,本底代表無模式的地方,飽和代表測量的數(shù)據(jù)超量程的地方,干擾代表噪聲.圖7 為該模式的頻譜圖,可以清晰地看到,這種模式的頻率范圍在80—120 kHz 之間.

圖6 聚類識別結(jié)果以及模式實際觀測到的位置 (a) 聚類識別結(jié)果;(b) 模式實際出現(xiàn)的位置;(c) 模式在托卡馬克中的位置Fig.6.Cluster recognition results and the position where the pattern is actually observed:(a) Cluster recognition results;(b) the position where the pattern actually appears;(c) the position of the pattern in the Tokamak.

圖7 對于64960 炮,模式頻率特征Fig.7.Mode frequency characteristics for shot 64960.

利用64960 炮識別出的該模式的典型信號進行模式篩選,在64962,64964,64965,64966,64967,64968,64969 炮也同樣發(fā)現(xiàn)了該模式.可以發(fā)現(xiàn)該模式存在于ECEI 第5 列附近,空間分布有一定特點;在頻譜圖上也具有一定特點.目前,還沒有對此類模式的記載,由此可以推斷模式很可能是一種新模式.值得注意的是,對于該模式的判定,以及它是否為新模式,仍需進一步物理上的分析.但譜聚類方法給出了一種尋找潛在的新模式的新思路,這在模式識別上具有很高的應(yīng)用價值.

4 結(jié)論

本文利用譜聚類的方法對EAST 裝置上不同診斷系統(tǒng)的數(shù)據(jù),包括ECEI、ECE、磁探針、SXR以及快輻射數(shù)據(jù)進行了分析,在識別精度以及效率方面表現(xiàn)良好.尤其可以對特定模式進行篩選,具有較大的實用性;此外填補了譜聚類算法在單點一維數(shù)據(jù)識別上的空白.由于在識別不同數(shù)據(jù)及不同模式時,算法本身不需要進行調(diào)整,因此表明其優(yōu)異的遷移性,為實際工作帶來了便利.利用譜聚類算法能為尋找潛在的新模式提供新思路,對等離子體物理的研究有很大的意義.目前定義數(shù)據(jù)之間的相似度使用的是距離度量方式,為進一步提高識別精度以及效率,下一步將會尋找更適合核聚變裝置數(shù)據(jù)的相似度度量方式.