基于彈道虧損形狀特征的堆積脈沖識別方法

王 明,周建斌,*,王懷平,2,汪雪元,2,劉 易,洪 旭

(1.成都理工大學 核技術與自動化學院,四川 成都 610059;2.東華理工大學,江西 南昌 330013)

X熒光光譜分析是一種速度快、非破壞性的測量方法,廣泛應用于現場礦物篩查分類、安檢等場景。X熒光光譜分析儀工作時,獲得穩定的光譜需要一定的測量時間,且元素含量越低,所需的測量時間越長,工程應用中一般通過提高X射線管的照射強度減少探測時間。核輻射探測系統中,脈沖信號衰減至基線需要一定的時間,在脈沖持續時間內,若其他輻射事件繼續與探測器響應,則會出現兩個或多個脈沖堆積的現象,形成畸變脈沖,也稱為堆積脈沖。堆積脈沖會導致能譜出現假峰(和峰)、本底提高、能量分辨率衰退等問題,常用數字脈沖成形技術及反堆積技術改善能譜性能[1]。高斯成形系統[2-3]及CR-RCn準高斯成形系統[4-7]是無限沖激響應系統,一般用于提升低計數率條件下能譜的能量分辨率。數字梯形成形算法能免疫彈道虧損、降低脈沖堆積的概率及提高信噪比[8-9],廣泛用于獲取高計數率與高能量分辨率能譜。此外,極零相消技術也是一種降低脈沖堆積概率的方法[10-11]。

快速成形技術是識別堆積脈沖的一種方法[12-15],該方法利用反卷積算法,將緩慢衰減的脈沖轉變成沖激脈沖。沖激脈沖時間短、脈沖通過率高,利用沖激脈沖的時間間隔可判斷梯形脈沖是否堆積,但該方法無法甄別時間間隔低于沖激脈沖分辨時間的堆積脈沖。此外,Imperiale等[16]利用梯形成形的上升時間來甄別堆積脈沖,Zhou等[17]采用梯形脈沖的寬度來甄別堆積脈沖,然而,梯形脈沖的上升時間及寬度隨著脈沖幅度增加而增加。Yu等[18]利用雙指數模型得到了平頂寬度穩定的梯形成形脈沖,并利用梯形平頂寬度來判別堆積脈沖,但測量時間特征受到ADC采樣精度的限制。Jordanov[19-20]提出了一種新的脈沖形狀甄別技術,該技術將原始脈沖成形為逆鋸齒波脈沖和梯形脈沖,計算出其幅度比值,該比值稱為時不變脈沖形狀標簽(TIPS),通過TIPS的大小判斷脈沖是否堆積。TIPS不受脈沖幅度變化及ADC分辨時間的影響,僅與彈道虧損相關,但逆鋸齒波成形器是高通濾波器,信噪比低,TIPS展寬大,不利于堆積脈沖識別。三角波成型器是帶通濾波器,信噪比高,TIPS展寬小,為此,本文提出梯形-三角脈沖形狀甄別技術,用于提升X熒光光譜測量系統檢測堆積脈沖的能力。

1 TIPS特征描述

利用數字成形技術將原始數字脈沖分別成形為梯形脈沖及三角脈沖。輻射事件響應的梯形脈沖及三角脈沖示意圖如圖1所示,其中,Ⅰ為非堆積,Ⅱ為兩個事件堆積,a為輻射事件序列,b、c分別為事件響應的梯形脈沖、三角脈沖。對于非堆積脈沖,梯形脈沖高度及三角脈沖高度均正比于輻射事件能量,TIPS是一個常數,與輻射事件能量無關。當脈沖堆積時,梯形脈沖能免疫彈道虧損,脈沖幅度依然正比于兩個輻射事件的能量。三角脈沖不能免疫彈道虧損,脈沖幅度損失,堆積脈沖的形狀特征TIPS變大。堆積時間間隔越大,TIPS越大。根據TIPS的大小,即可判斷脈沖是否發生堆積。

圖1 梯形脈沖及三角脈沖成形Fig.1 Trapezoidal pulse and triangular pulse

(1)

其中,PKflat、PKpin分別為梯形脈沖高度及三角脈沖高度。

2 原理

2.1 TIPS特征分析

經前端電路處理后的核脈沖信號,首先快速上升至峰值,然后緩慢衰減至基線,常用雙指數衰減脈沖信號x(t)進行模擬:

(2)

其中:A為比例常數,與探頭和電子學電路相關;E為輻射事件的能量;u(t)為單位階躍信號,t為時間;τ及θ為衰減常量,τ與C-R微分電路相關,θ與探頭和電路特性相關,θ通常遠大于τ。τ描述的緩慢衰減指數信號用于表征核脈沖的慢成份;θ描述的快速衰減指數信號用于表征核脈沖的快成份。

梯形(三角)脈沖信號為:

(t-tb)u(t-tb)+(t-tc)u(t-tc)]

(3)

其中:yT為梯形脈沖信號;bT為比例常數;ta為上升時間,tb-ta為平頂脈沖信號時間,tc=ta+tb。當t在[ta,tb]之間時,最大值為bTE。當ta=tb時,式(3)可描述為等腰三角脈沖信號,ya為三角脈沖信號,當t=ta時,最大值為baE,ba為比例常數。

對于兩個輻射事件堆積的情況,假設首先與探測器響應的輻射事件能量為El,間隔時間為ΔT時,能量為Er的輻射事件與探測器響應,且ΔT

ya(t) =ElV(t) +ErV(t-ΔT)

(4)

其中:tw為三角脈沖的寬度;V為三角脈沖信號;ya為重疊的三角脈沖信號。

設等腰三角波的寬度為tw,即tw=tc=2ta=2tb,堆積三角波的不可導位置t=tw/2或t=ΔT+tw/2的極值分別為ba(El+Er-2ΔTEr/tw)、ba(El+Er-2ΔTEl/tw)。當ba(El+Er-2ΔTEr/tw)≥ba(El+Er-2ΔTEl/tw),即El≥Er時,極大值點為ba(El+Er-2ΔTEr/tw);當ba(El+Er-2ΔTEr/tw)ΔT,堆積的梯形脈沖幅度為bT(El+Er)。

根據式(1),非堆積脈沖的TIPS為:

(5)

其中,TIPSa為非堆積脈沖的梯形-三角波形狀特征,簡稱為三角波形狀特征。

兩個事件堆積時的TIPS為:

(6)

圖2 TIPS-El/Er曲線Fig.2 TIPS-El/Er curve

2.2 甄別器評估方法

脈沖形狀特征甄別器是一種二分類算法,即對非堆積脈沖和堆積脈沖進行分類。常用精確率、召回率以及F1得分來度量分類器的性能。精確率與召回率是矛盾的,精確率越高,召回率就越低。F1得分是精確率與召回率的調和平均數,它同時兼顧了分類模型的精確率與召回率,計算公式分別為:

(7)

(8)

(9)

其中:Pacc為精確率;Pcal為召回率;TP為正確的正例;FN為錯誤的反例;FP為錯誤的正例。本文將堆積脈沖設為正例,非堆積脈沖設為反例。精確率表示甄別系統正確識別堆積脈沖的比例;召回率表征甄別系統識別到的堆積脈沖與全部堆積脈沖的比值。

3 數字成形算法實現

對式(2)進行ADC采樣量化,獲得數字脈沖信號,采樣時間為Ts。對數字脈沖進行z變換得z變換表達式為:

(10)

梯形(三角脈沖)傳遞函數HT(z)[17]為:

(11)

其中:na=ta/Ts;nb=tb/Ts;nc=tc/Ts。當na=nb時可得到三角波成形傳遞函數Ha(z)。

在硬件上實現數字成形算法以整數運算為基礎,不存在溢出的情況下進行加法、減法以及乘法運算不會產生計算誤差,除法運算容易產生截斷誤差。采用級聯結構來設計數字成形算法可降低系統的復雜度。對于線性時不變系統,交換級聯的順序不影響最終的結果輸出。差分項放在前級,累加項放在后級,能夠防止溢出。

為計算TIPS,需根據式(11)將原始脈沖成形為三角脈沖和梯形脈沖。通過反卷積算法將原始脈沖x成形為沖激脈沖v,如圖3a所示。為了避免浮點運算,將d1和d2改為有理分式形式。設d1≈nu1/de1、d2≈nu2/de2。其中nu1、de1、nu2、de2均為整數。級聯項1-d1z-1及1-d2z-1分別改為(de1-nu1z-1)/de1以及(de2-nu2z-1)/de2。設三角脈沖的上升時間為na,縮放系數設為k1=D1/2M1,D1、M1均為整數。三角成形系統的硬件算法如圖3b所示,ya為三角脈沖信號。設梯形脈沖的上升時間為na1,平頂時間為nb1-na1,縮放系數設為k2=D2/2M2,D2、M2均為整數。梯形成形系統的硬件算法如圖3c所示,yT為梯形脈沖信號。

圖3 算法實現邏輯圖Fig.3 Logic diagram of algorithm implementation

算法部署于xc7z020CLG400-1芯片中,參數取值分別為:de1=100,nu1=95,de2=100,nu2=61,na=4,na1=8,nb1=12,D1=1,D2=800,M1=16,M2=17。算法資源消耗為:19個LUT,7個LUTRAM及359個FF。

4 實驗

利用基于FAST-SDD探測器的CIT-3000MD X射線熒光光譜測量系統進行實驗。FAST-SDD探測器(XR-100SDD)由Amptek生產。XR-100SDD內部封裝了熱電冷卻固態探測器和前置放大器。XR-100SDD引出的獨立電源線和信號線接入新先達(M&C Co. Ltd)的數字脈沖處理系統。光源使用科頤維公司設計的Ag靶X光管(KYW2000A),額定工作電壓與額定工作電流分別為50 kV、1 mA。

梯形成形脈沖的上升時間為400 ns,平頂寬度為200 ns,脈沖持續時間為1 μs。三角脈沖的寬度為350 ns,沖激脈沖的時間分辨率約為150 ns(20 MHz 采樣率條件下)[14]。采用鉛黃銅樣品進行實驗,測量時間為100 s,光管的工作電壓為33 kV,工作電流分別為3.9 μA、1 mA,沖激脈沖通過率分別為6.5 kcps/s、1.22 Mcps/s。

4.1 TIPS分布圖

由于噪聲以及彈道虧損等的影響,非堆積脈沖的TIPS并非固定值,而是符合一定分布,通過TIPS分布圖可確定分割閾值。通過調整脈沖幅度的比例,使得TIPS譜線的均值為1 000,圖4a為管流為3.9 μA時的TIPS分布,Pi為逆鋸齒波形狀特征TIPS分布[20]、Pa為三角波形狀特征TIPS分布。低計數率時,TIPS譜線近似為非堆積脈沖的TIPS譜線,Pa展寬小,Pi展寬大。圖4b為管流為1 mA時的TIPS分布。高計數率時,由于脈沖堆積的影響,TIPS分布高值區域出現長的拖尾。在高值區長尾范圍內,隨TIPS的增加,Pi的概率密度快速降低,Pa的概率密度緩慢衰減。

圖4 TIPS分布Fig.4 TIPS distribution

表1 不同α對應的甄別閾值xαTable 1 Discrimination threshold xα corresponding to different α

4.2 測量結果及討論

SPT譜中,銅元素Kα、Kβ特征X射線能量位于(7.6 keV,9.5 keV)能區。當脈沖通過率為6 500 s-1時,幾乎沒有堆積脈沖,高能區平坦。Kα峰與Kβ峰面積占總面積的95.71%。兩個脈沖堆積時,一階和峰位于(15 keV,18 keV)能區;3個脈沖堆積時,二階和峰位于(20 keV,26 keV)能區;3個以上脈沖堆積的概率較低,且已超出ADC的測量范圍,可忽略。

管流為1 mA時,梯形成形能譜如圖5所示。SPT是利用傳統的快成形技術[14]甄別堆積梯形脈沖獲得的能譜。在快成形技術甄別堆積脈沖的基礎上,利用梯形-逆鋸齒波形狀特征[20]或梯形-三角波形狀特征進一步甄別堆積脈沖,分別獲得能譜ISA-SPT[20]及TRI-SPT,其中非堆積脈沖損失率為15%。SPT中一階和峰、二階和峰計數率分別為26 580、2 030 s-1,ISA-SPT中一階和峰、二階和峰計數率分別為8 400、181 s-1,TRI-SPT中一階和峰、二階和峰計數率分別為6 750、116 s-1。因此,使用梯形-三角波形狀特征甄別堆積脈沖,可進一步降低和峰計數率。

圖5 拒絕堆積脈沖的能譜對比Fig.5 Comparison of energy spectrum rejected pile-up pulse

取Cu元素的Kα、Kβ峰射線全能峰以及其和峰計算出的精確率、召回率、F1得分列于表2。由表2可知:當α為5%時,F1得分最高;隨α的增大,召回率逐漸增大,精確度逐漸降低;相同α條件下,三角波TIPS甄別法識別率高于逆鋸齒波TIPS甄別法[20]。當α為15%時,三角波TIPS甄別法識別的精確率、召回率、F1得分分別為73.55%、78.75%、76.06%,分別提升了6.39%、7.66%、7.08%。

表2 不同α對應的堆積脈沖分類結果Table 2 Classification result of pile-up pulse corresponding to different α

峰總比表示全能峰的探測效率,可評價能譜的質量,峰總比表示全能峰面積與總面積的比。當電流為3.9 μA時,Cu元素Kα峰的峰總比為82.13%;當光管電流為1 mA時,Kα峰的峰總比降到62.93%。利用形狀特征TIPS識別并拒絕堆積脈沖后,峰總比的測量結果如圖6所示。PTi、PTa分別表示ISA-SPT、TRI-SPT能譜中Cu元素Kα峰的峰總比,峰總比隨α的增加而增加。當α相同時,PTa高于PTi。當α=15%時,PTi、PTa分別為74.86%、76.60%,峰總比提升了1.74%。

圖6 峰總比的測量結果Fig.6 Measurement result of peak to total ratio

5 結論

本文提出了一種基于彈道虧損形狀特征的堆積脈沖甄別技術,其原理為:首先將原始脈沖成形為三角脈沖和梯形脈沖;其次計算其幅度比值,即脈沖形狀特征值TIPS;最后根據TIPS判斷脈沖是否發生堆積。由于TIPS與脈沖幅度及ADC采樣率無關,在20 MHz采樣率條件下,仍具有比較好的堆積脈沖甄別能力。通過理論計算分析了兩個脈沖堆積時的TIPS特征值,堆積時間間隔越大或堆積事件的能量越接近,TIPS越大,堆積脈沖越容易甄別。

使用提出的方法識別鉛黃銅樣品中1.22×106s-1脈沖通過率的堆積脈沖。當非堆積脈沖的損失率為15%時,識別堆積脈沖的精確率、召回率以及F1得分分別為73.55%、78.75%以及76.06%,與文獻[20]提出的TIPS甄別技術進行對比,分別提升了6.39%、7.66%以及7.08%。Cu元素的Kα峰的峰總比為76.60%,提升了1.74%。該方法減少了偽峰的計數率,在現場礦物篩查分類以及安檢等快速識別微量元素場景中具有潛在優勢。