趙子博 莊革 謝錦林 渠承明 強(qiáng)子薇

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

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

1 引言

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

核聚變裝置已經(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)電子回旋輻射計(jì)(ECE)[4]和 ECEI[5?9]等,因此總數(shù)據(jù)量巨大.然而目前主要依賴人工進(jìn)行模式區(qū)分,該方法效率低、成本高,無法滿足實(shí)際需求,導(dǎo)致大量的診斷數(shù)據(jù)沒有得到充分分析.因此尋找高效模式識別的方法十分重要.

近年來,機(jī)器學(xué)習(xí)在核聚變領(lǐng)域已經(jīng)有了廣泛而深入的應(yīng)用.機(jī)器學(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í)需要大量帶有標(biāo)簽的數(shù)據(jù)來進(jìn)行訓(xùn)練,而目前大多數(shù)診斷數(shù)據(jù)的標(biāo)簽尚未完善;另一方面監(jiān)督學(xué)習(xí)的遷移性很差,對于不同裝置或者同一裝置上的不同診斷系統(tǒng),甚至同一裝置上的同一個(gè)診斷系統(tǒng)不同條件下的數(shù)據(jù)都需要重新訓(xùn)練模型,效率低、適用性差.而無監(jiān)督學(xué)習(xí)不存在此缺點(diǎn),因此本文采用無監(jiān)督學(xué)習(xí)的方法.

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

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

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

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

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

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

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

進(jìn)一步,定義矩陣:

定義

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

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

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

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

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

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

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

圖3 各信道的信號時(shí)序圖 (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)定性演化過程圖,標(biāo)號(1)—(8) 依次對應(yīng)8 個(gè)時(shí)刻點(diǎn);圖4(b)黑色、紅色、藍(lán)色曲線分別代表圖4(a)各圖對應(yīng)顏色點(diǎn)處的時(shí)序圖.從圖4(b)可以明顯看出:最開始鋸齒崩塌結(jié)束,冷磁島占據(jù)整個(gè)q=1 面;隨著等離子體加熱,鋸齒爬升,芯部電子溫度緩慢提高,之后重聯(lián)發(fā)生,芯部熱量向外輸運(yùn);最后鋸齒崩塌,冷磁島重新占據(jù)整個(gè)q=1 面,符合鋸齒模的演化過程.從圖4(a)可以看出,整個(gè)演化圖的空間結(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)隨時(shí)間的演化過程;(b) 黑色、紅色、藍(lán)色曲線分別代表圖4(a)各圖對應(yīng)顏色點(diǎn)處的時(shí)序圖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).

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

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

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

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

對ECE 診斷的時(shí)序信號進(jìn)行傅里葉變換,得到各個(gè)時(shí)間點(diǎn)的頻率信息,每個(gè)時(shí)間點(diǎn)的頻率序列對應(yīng)一個(gè)樣本,同種模式的表現(xiàn)特征是頻率序列強(qiáng)度相似.各個(gè)樣本的相似度用樣本之間的指數(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)有展寬非常大的破裂,符合鋸齒模的特征,認(rèn)定識別出的模式為鋸齒模.對于SXR、快輻射以及ECE 三種診斷數(shù)據(jù)仿照前述操作進(jìn)行時(shí)間聚類,對50015—50115 炮的24000 個(gè)時(shí)間片段的識別結(jié)果進(jìn)行統(tǒng)計(jì),其中6700 個(gè)時(shí)間片段被識別為鋸齒模.根據(jù)頻譜圖,可以判定有6555 個(gè)真實(shí)存在鋸齒模的時(shí)間片段.同時(shí),根據(jù)頻譜圖也可以判定出在被聚類算法識別出的6700 個(gè)鋸齒模片段中,有6057 個(gè)是真實(shí)的鋸齒模片段.因此,根據(jù)(20)式和(21)式對查準(zhǔn)率和查全率的定義,可以計(jì)算得到P=90.4%,R=92.4%.

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

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

3.3 新模式探索

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

圖6 聚類識別結(jié)果以及模式實(shí)際觀測到的位置 (a) 聚類識別結(jié)果;(b) 模式實(shí)際出現(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 炮識別出的該模式的典型信號進(jìn)行模式篩選,在64962,64964,64965,64966,64967,64968,64969 炮也同樣發(fā)現(xiàn)了該模式.可以發(fā)現(xiàn)該模式存在于ECEI 第5 列附近,空間分布有一定特點(diǎn);在頻譜圖上也具有一定特點(diǎn).目前,還沒有對此類模式的記載,由此可以推斷模式很可能是一種新模式.值得注意的是,對于該模式的判定,以及它是否為新模式,仍需進(jìn)一步物理上的分析.但譜聚類方法給出了一種尋找潛在的新模式的新思路,這在模式識別上具有很高的應(yīng)用價(jià)值.

4 結(jié)論

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