基于改進爬山法的帶電粒子圓徑跡重建

牛莉博，李玉蘭，傅楗強，黃 孟，何 彬，李元景

（1.清華大學 工程物理系，北京 100084；2.第二炮兵工程大學102教研室，陜西 西安 710025；3.清華大學 粒子技術(shù)與輻射成像教育部重點實驗室，北京 100084）

徑跡探測器是一重要的電離輻射探測器，被廣泛應用于大型物理實驗、核技術(shù)及應用等領(lǐng)域，基于時間投影室（TPC）的徑跡探測器是目前使用最為廣泛的徑跡探測器之一［1-3］。TPC的基本原理為：當帶電粒子經(jīng)過TPC 靈敏體積時，使TPC 工作氣體電離產(chǎn)生電子-離子對，電子在TPC 內(nèi)部均勻電場作用下，向讀出端蓋漂移并產(chǎn)生雪崩放大。電子信號由按一定規(guī)則排列的若干讀出盤讀出，讀出盤的位置信息確定了帶電粒子x-y 方向的位置，而z 方向的位置由z＝vD（t1-t2）確定（vD為電子在電場中的漂移速度；t1為帶電粒子穿過TPC 的時間，由外部時間探測器的觸發(fā)信號給出；t2為電子到達讀出端面的時間）。根據(jù)讀出盤的擊中信息重建帶電粒子的飛行徑跡是TPC 數(shù)據(jù)分析的重要步驟，在均勻磁場情況下，帶電粒子的徑跡是1條空間螺旋線，其在二維讀出平面上的投影為圓形。圓徑跡重建的任務為：對于給定的1 組二維讀出平面上的擊中點（xi，yi）（i＝1，2，…，n），通過一定的算法判斷這些擊中點是否來自于同一圓徑跡，即這些擊中點是否由同一帶電粒子產(chǎn)生；若這些測量點由同一帶電粒子產(chǎn)生，則對這些擊中點進行擬合得到該圓徑跡參數(shù)，即圓心（a，b）和半徑R 的擬合值。

霍夫變化是一種簡單有效的徑跡重建方法，在粒子徑跡重建算法中得到了廣泛使用［4-6］。然而，由于霍夫變換需將參數(shù)空間進行離散化，隨參數(shù)空間維數(shù)的增加，計算量將急劇增加，計算速度變慢。文獻［7］指出，霍夫變換法已很難滿足對圓徑跡的識別與檢測。爬山法最早由Yager和Filev［8-9］提出，是一有效的快速聚類算法。與霍夫變化一樣，傳統(tǒng)的爬山法也需將參數(shù)空間進行網(wǎng)格化，對高維數(shù)據(jù)的聚類速度急劇下降。而改進爬山法則使用數(shù)據(jù)點替代網(wǎng)格點，在數(shù)據(jù)點構(gòu)造密度函數(shù)替代在網(wǎng)格點構(gòu)造山峰函數(shù)，通過尋找密度函數(shù)最大值以確定參數(shù)擬合結(jié)果，可極大提高計算速度，對高維數(shù)據(jù)的聚類分析具有良好適應性。本文通過空間映射將圓徑跡尋跡問題轉(zhuǎn)化為尋找參數(shù)空間上的聚類中心問題，使用改進爬山法獲得圓徑跡的參數(shù)擬合值，并給出不同情況下圓徑跡擬合參數(shù)的精度。

1 方法

1.1 爬山法

假設有n個數(shù)據(jù)點｛x1，x2，…，xn｝，每個數(shù)據(jù)點是p 維的，即xj∈Rp（j＝1，…，n）。傳統(tǒng)爬山法中，首先在數(shù)據(jù)空間中定義1組網(wǎng)格點｛Ni｝，Ni∈RF，聚類中心被固定在這些網(wǎng)格點上。然后在每個網(wǎng)格點上定義山峰函數(shù)如下：

其中，d（xj，Ni）為數(shù)據(jù)點xj和網(wǎng)格點Ni之間的距離。顯然，對于擁有較多臨近點的網(wǎng)格點，其山峰函數(shù)值較大，參數(shù)α決定臨近區(qū)域的大小。具有大的山峰函數(shù)值M1（Ni）的網(wǎng)格點Ni成為聚類中心的潛在可能性大。于是，可將具有最大峰值的網(wǎng)格點設定為第1個聚類中心，即可通過以下目標函數(shù)找到第1個聚類中心N＊1：

對于具有多個聚類中心的情況，在尋找下一聚類中心前，有必要消除上一網(wǎng)格點對下一聚類中心的影響。在找到第k-1個聚類中心N＊k-1后，估計第k個聚類中心的計算公式如下：

其中，k＝2，3，…，且有：

式（4）為峰值調(diào)整方程，參數(shù)γ決定了消減已提取聚類中心點峰值影響的臨近區(qū)域的大小。通過式（3）和（4）的迭代可得到新的聚類中心。

1.2 改進爬山法

顯然，爬山法的計算時間隨數(shù)據(jù)維度的增加而急劇增加，Chui［10］用數(shù)據(jù)點代替網(wǎng)格點來改進原始的爬山算法。在此基礎上，Yang等［11］對山峰函數(shù)和峰值調(diào)整方程進行了修改，進一步減小算法的計算時間。在改進的爬山算法中，考慮用核密度估計每個數(shù)據(jù)點向xi的密度：

且

式中：參數(shù)ρ為歸一化項；參數(shù)λ可由自相關(guān)有關(guān)算法得到，或由實驗者給定。于是第1個聚類中心估計公式為：

在估計出第k-1個聚類中心后，用改進的山峰調(diào)整算法調(diào)整下一輪的山峰：

在調(diào)整后的密度中，具有最高密度值Pk（xi）的那個數(shù)據(jù)點即為所要提取的第k個聚類中心：

通過式（8）和（9）的不斷迭代，可提取出新的聚類中心。

1.3 空間映射

盡管改進爬山法能很好估計出聚類中心，但由于初始測量點并不具備聚類特性，無法直接應用于圓徑跡的重建。為能將改進爬山法應用于圓徑跡的重建，首先需對擊中點進行變換，即將測量得到的二維坐標數(shù)據(jù)（xi，yi）轉(zhuǎn)變?yōu)閰?shù)空間中的參數(shù)值，該參數(shù)值與對應圓徑跡參數(shù)有關(guān)。圓方程可如下表示：

式中，d＝R2-a2-b2。由于平面上任意3個不共線的點將確定1個圓，且該圓的圓心和半徑可按下式計算得到：

這樣，對所有（xi，yi）中任意3個不共線的測量點組合便可得到個參數(shù)值（aj，bj，Rj）。由于測量點（xi，yi）來自于相同的圓徑跡，因此由測量點空間映射得到的參數(shù)值（aj，bj，Rj）在參數(shù)空間里面形成聚類。圖1a為有兩條圓徑跡并包含噪聲點的原始測量數(shù)據(jù)分布，圖1b、c、d分別為測量點空間映射后所得參數(shù)值在a-b、a-R、b-R 平面上的投影。由圖1可知，參數(shù)值在參數(shù)空間中存在聚類中心，且該聚類中心在參數(shù)空間中的位置即為圓徑跡參數(shù)的擬合值。

2 算法性能

2.1 單圓徑跡重建

計算中選取的圓徑跡的真實參數(shù)為圓心（a，b）＝（4，5），半徑R＝2，使用該參數(shù)隨機產(chǎn)生的擊中點數(shù)目分別為8、12和16。在實際探測器中，測量點（xi，yi）的坐標存在測量誤差，模擬中隨機產(chǎn)生的測量點（xi，yi）到圓心的距離服從高斯分布，該分布的均值為半徑R，方差分別為R的5%、10%和15%。除此之外，模擬數(shù)據(jù)還隨機產(chǎn)生了一些噪聲點，噪聲點數(shù)分別選為數(shù)據(jù)點的50%、25%和0%，噪聲點的（xi，yi）坐標均勻分布在［0，10］區(qū)間上。圖2為上述模擬數(shù)據(jù)中的1個事例，圖2a為隨機產(chǎn)生的原始數(shù)據(jù)，有效數(shù)據(jù)點為20個，噪聲數(shù)據(jù)點為10個，即50%的噪聲數(shù)據(jù)點，（xi，yi）坐標的測量不確定度為10%的R；圖2b為數(shù)據(jù)點空間映射產(chǎn)生的參數(shù)值在參數(shù)空間a-b平面上的投影；圖2c為密度函數(shù)值在參數(shù)空間a-b 上的投影；圖2d為參數(shù)空間上的密度函數(shù)值，密度函數(shù)值的大小用彩虹圖表示，數(shù)值由大到小依次由紅色到藍色表示。

式（5）中參數(shù)λ 的選擇決定了改進爬山法是否能正確獲得擬合參數(shù)以及參數(shù)的擬合精度，文獻［11］指出λ 的選擇可通過自相關(guān)的有關(guān)算法給出或由實驗者設定。實際上參數(shù)λ的選擇與參數(shù)空間中的數(shù)據(jù)結(jié)構(gòu)有直接關(guān)系，模擬研究發(fā)現(xiàn)對于上述模擬數(shù)據(jù)，當λρ＝50時，改進爬山法具有良好的適用性，因此本文所有計算均選擇λρ＝50。在實際的徑跡重建中，應當根據(jù)探測器的尺寸、圓徑跡的分布及噪聲水平等探測器性能參數(shù)合理設置λ。計算中，在每種情況下隨機產(chǎn)生10 000條模擬徑跡，使用改進爬山法得到每條徑跡的擬合參數(shù)，表1列出了在不同情況下改進爬山法給出的圓徑跡參數(shù)的擬合結(jié)果及其精度。由計算結(jié)果易知，在不同的噪聲水平和測量精度下，改進爬山法均能給出較高精度的圓徑跡參數(shù)擬合值。但隨測量不精確度和噪聲水平的增加，參數(shù)的擬合精度在下降。

圖1 測量點及其空間映射Fig.1 Original data points and their transformations in parameters'space

圖2 改進爬山法對單圓徑跡擬合結(jié)果Fig.2 Fitting results of single circle track based on modified mountain method

表1 不同情況下改進爬山法對圓徑跡參數(shù)擬合精度Table 1 Fitting values and their errors derived by modified mountain method under various scenarios

2.2 其他圓徑跡重建

圖3為高能物理實驗中帶電粒子在TPC中的投影徑跡，TPC 的中心為中心束流管。高能粒子在束流管中的對撞中心點發(fā)生碰撞，產(chǎn)生的初級帶電粒子的徑跡在讀出平面上的投影為1段圓弧，該圓弧通過坐標原點，其所對應的圓心角由粒子的動量和磁場強度決定；而對于次級帶電粒子的徑跡，其圓弧型徑跡并不一定通過坐標原點，對應的圓心角同樣由粒子的動量和磁場強度決定。

圖3 帶電粒子在TPC中的投影徑跡Fig.3 Track projection of particles'trajectories in TPC

針對幾種不同的圓弧型徑跡，圖4示出了使用改進爬山法的重建結(jié)果。圖4a～c為單圓弧型徑跡的重建結(jié)果，圖4a為過原點的初級帶電粒子產(chǎn)生的圓弧型徑跡，該徑跡完全投影到了讀出平面，圖4b為僅部分徑跡投影到TPC讀出平面的情況，而圖4c則為次級粒子產(chǎn)生的徑跡情況。圖4d～f為上述幾種情況下的雙圓徑跡組合情況下的徑跡重建結(jié)果。與單圓徑跡重建相同，計算中使用參數(shù)λρ＝50，噪聲點數(shù)取信號點數(shù)的50%～100%，位置測量的隨機誤差為徑跡半徑R 的5%，每條徑跡上信號點數(shù)在16～21之間。計算結(jié)果顯示，改進爬山法同樣適用于多條圓徑跡的重建。

3 結(jié)語

本文通過空間映射將二維平面上的帶電粒子圓徑跡重建問題轉(zhuǎn)變?yōu)閰?shù)空間上的聚類問題，其中參數(shù)聚類中心在參數(shù)空間上的位置代表了圓徑跡的參數(shù)，使用改進爬山法在參數(shù)空間上尋找該聚類中心。充分考慮探測器的位置測量不精確度、噪聲水平和測量點數(shù)等因素，在不同情況下使用改進爬山法計算得到了圓徑跡參數(shù)的擬合值及其精度。計算結(jié)果顯示，該方法具有較強的魯棒性和抗噪聲干擾能力，在不同噪聲水平下均能給出較高的擬合精度。此外，使用改進爬山法對二維平面上的兩條圓徑跡進行了重建，結(jié)果顯示該算法同樣適用多條圓徑跡的重建。

圖4 不同類型的圓徑跡重建Fig.4 Track reconstruction of circles under various scenarios

［1］ HILKE H J.Time projection chambers［J］.Reports on Progress in Physics，2010，73（11）：116201-116236.

［2］ BOWDEN N S，HEFFNER M，CAROSI G，et al.Directional fast neutron detection using a time projection chamber［J］.Nuclear Instruments and Methods in Physics Research A，2010，624（1）：153-161.

［3］ HANG M，LI Y L，DENG Z，et al.A fast neutron spectrometer based on GEM-TPC［C］∥IEEE Nuclear Science Symposium and Medical Imaging Conference Record （NSS／MIC）.［S.l.］：［s.n.］，2012：146-148.

［4］ FABBIETTI L，ANGERER H，ARORA R，et al.The PANDA GEM-based TPC prototype［J］.Nuclear Instruments and Methods in Physics Research A，2011，628（1）：204-208.

［5］ CHESHKOV C.Fast Hough-transform track reconstruction for the ALICE TPC［J］.Nuclear Instruments and Methods in Physics Research A，2006，566（1）：35-39.

［6］ STRANDLIE A，F(xiàn)RUHWIRTH R.Track and vertex reconstruction：From classical to adaptive methods［J］.Reviews of Modern Physics，2010，82（2）：1 419-1 458.

［7］ TSUJI S，MATSUMOTO F.Detection of ellipses by a modified Hough transformation［J］.IEEE Transactions on Computers，1978，27（8）：777-781.

［8］ YAGER R R，F(xiàn)ILEV D P.Approximate clustering via the mountain method［J］.IEEE Transactions on System，Man and Cybernetics，1994，24（8）：1 279-1 284

［9］ YAGER R R，F(xiàn)ILEV D P.Generation of fuzzy rules by mountain clustering［J］.Journal of Intelligent and Fuzzy Systems，1994，2（3）：209-219.

［10］CHUI S L.Fuzzy model identification based on cluster estimation［J］.Journal of Intelligent and Fuzzy Systems，1994，2（3）：267-278.

［11］YANG M S，WU K L.A modified mountain clustering algorithm［J］.Pattern Analysis and Applications，2005，8（1）：125-138.