唐 琳,趙衛(wèi)東,余松科,劉 澤,余小東,孟 源,黃興祿
1. 成都大學(xué)電子信息與電氣工程學(xué)院,四川 成都 610106 2. 數(shù)學(xué)地質(zhì)四川省重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué)),四川 成都 610059 3. 成都大學(xué)科研處,四川 成都 610106 4. 成都大學(xué)計(jì)算機(jī)學(xué)院,四川 成都 610106
輻射測(cè)量中最重要的兩個(gè)指標(biāo),一個(gè)是能量分辨率,另一個(gè)則是計(jì)數(shù)率。但能量分辨率和計(jì)數(shù)率是相互矛盾的,改善能量分辨率將降低計(jì)數(shù)率,相反提高計(jì)數(shù)率將損失能量分辨率[1-3]。探測(cè)器的能量分辨率由多個(gè)因素決定,比如探測(cè)信號(hào)的統(tǒng)計(jì)漲落、信號(hào)處理器的噪聲、外部干擾以及溫度等。統(tǒng)計(jì)漲落決定了給定探測(cè)器能量分辨率的理論極限,而其他因素的影響則可以通過適當(dāng)?shù)脑肼暈V除和電子技術(shù)來降低[4-6]。
以往關(guān)于能量分辨率的研究大多利用譜反卷積對(duì)獲取到的能譜進(jìn)行后處理,從而降低特征峰的半高寬。這些后處理方法是通?;趯@取到的能譜建模為輸入能譜和探測(cè)器響應(yīng)函數(shù)這兩個(gè)隨機(jī)變量的函數(shù),而探測(cè)器的響應(yīng)函數(shù)又可以建模為輸入脈沖高度與輸出脈沖高度的聯(lián)合概率分布。目前至少存在三種對(duì)譜線進(jìn)行后處理的方法可以提高能量分辨率,包括使用正則化、極大似然和最大熵。然而,這三種后處理方法都涉及了為光譜反褶積建模探測(cè)器的能量響應(yīng)函數(shù),因此這些方法將是一項(xiàng)計(jì)算量非常大的任務(wù)[7]。
傳統(tǒng)的X射線光譜數(shù)據(jù)處理平臺(tái)采用的是多道譜的成譜方法,其主要特征在于譜圖中每一個(gè)道址上計(jì)數(shù)的增量是基于單個(gè)脈沖幅度的,也就是說,每一個(gè)脈沖幅度被接收,相應(yīng)的道址上的計(jì)數(shù)值就會(huì)加一。早在2005年,國外就有研究人員初步提出了多脈沖成譜的相關(guān)方法[8],但國內(nèi)尚無相關(guān)成果公布。筆者在前期研究工作[9-10]中已經(jīng)提出了脈沖剔除法和脈沖修復(fù)法處理突變脈沖得到真實(shí)可靠的脈沖幅度信息,并因此降低了統(tǒng)計(jì)漲落。本文在上述研究基礎(chǔ)上提出一種X射線光譜數(shù)據(jù)處理平臺(tái)的優(yōu)化設(shè)計(jì),采用可在線實(shí)施的多脈沖局部平均(MPLA)成譜技術(shù),選擇一個(gè)恰當(dāng)范圍的平均窗口對(duì)多個(gè)脈沖幅度進(jìn)行平均,平均后的幅度值被接收,再將相應(yīng)的道址上的計(jì)數(shù)值加一,該技術(shù)的實(shí)施使得X射線光譜數(shù)據(jù)處理平臺(tái)得以優(yōu)化。本文通過理論推導(dǎo)揭示了MPLA概率密度變換的基本理論特征,也通過實(shí)驗(yàn)驗(yàn)證了MPLA算法的特征,證明了在具有正態(tài)分布PDF的頻譜峰值的典型情況下,即使僅對(duì)兩個(gè)脈沖高度進(jìn)行平均,變換后峰的FWHM也會(huì)變窄。
MPLA算法與傳統(tǒng)多道譜的根本區(qū)別在于前者使用多個(gè)脈沖高度的平均值來確定需要增加計(jì)數(shù)的道址,而后者則是每讀取一個(gè)脈沖幅度值,就給對(duì)應(yīng)道址上的計(jì)數(shù)加一。MPLA算法的核心在于多個(gè)脈沖幅度的平均,但這里的平均并不是針對(duì)所有的脈沖幅度。如果對(duì)所有測(cè)量到的脈沖幅度都進(jìn)行平均,那么最終得到的譜圖就會(huì)完全失去其原有的統(tǒng)計(jì)特征,也會(huì)破壞不同元素特征峰之間的計(jì)數(shù)差異,因此MPLA算法中的平均必須是針對(duì)選定窗口內(nèi)一定脈沖幅度數(shù)量進(jìn)行平均,再通過平均得到的脈沖幅度值來確定增加計(jì)數(shù)的道址。
MPLA算法涉及兩項(xiàng)可變參數(shù),一是平均窗口的大小,本文以參數(shù)r表示,另一項(xiàng)參數(shù)則是每一次平均的脈沖幅度數(shù)量,本文以參數(shù)n表示。在下文對(duì)MPLA算法原理的描述中,我們?nèi)∑骄翱趓等于2,每一次參與平均的脈沖幅度數(shù)量為2。
在MPLA算法的執(zhí)行流程包含以下幾個(gè)步驟,如圖1所示。首先讀取第一個(gè)脈沖幅度時(shí)定位一個(gè)平均窗口,當(dāng)r為2時(shí)窗口內(nèi)包含的道址數(shù)為5,脈沖幅度讀取成功后更新當(dāng)前平均窗口的幅度和脈沖個(gè)數(shù),每次更新后即對(duì)平均窗口內(nèi)的脈沖個(gè)數(shù)進(jìn)行判斷,當(dāng)其小于預(yù)設(shè)的參數(shù)n時(shí),則繼續(xù)讀取下一個(gè)脈沖幅度,當(dāng)平均窗口內(nèi)的脈沖幅度數(shù)量等于參數(shù)N時(shí)則對(duì)相應(yīng)平均窗口內(nèi)的脈沖幅度進(jìn)行平均,得出的平均數(shù)即為需要更新計(jì)數(shù)的道址,然后再對(duì)取平均值的窗口內(nèi)脈沖幅度和脈沖個(gè)數(shù)進(jìn)行清零。
為了詳細(xì)闡述MPLA算法的實(shí)現(xiàn)過程,下文將以圖文結(jié)合的方式呈現(xiàn)一個(gè)MPLA實(shí)現(xiàn)示例。在當(dāng)前示例中,MPLA算法的平均窗口r設(shè)置為2,每一次參與平均的脈沖幅度數(shù)量設(shè)置為2。在圖2所示的MPLA算法步驟一中,讀取到的第一個(gè)脈沖幅度值為660,平均窗口的中心則為660,而窗口大小為2則代表窗口的范圍被限定在658~622之間。
此時(shí)活動(dòng)平均窗口僅有一個(gè),并且該窗口內(nèi)可用于平均的脈沖個(gè)數(shù)小于2,開始執(zhí)行MPLA算法示例的步驟二,如圖3所示,本次讀取的脈沖幅度為650,平均窗口中心移到650,平均窗口范圍為648~652,該范圍與前一個(gè)脈沖所處的范圍并無重疊,這就產(chǎn)生了一個(gè)新的平均窗口,對(duì)該窗口內(nèi)的脈沖幅度和以及脈沖個(gè)數(shù)進(jìn)行更新后再次判斷是否有達(dá)到平均條件的窗口,若沒有繼續(xù)執(zhí)行步驟三。
圖1 MPLA算法實(shí)現(xiàn)流程Fig.1 MPLA algorithm implementation flow
圖2 MPLA算法實(shí)現(xiàn)示例步驟一Fig.2 Step 1 of MPLA algorithm implementation example
在執(zhí)行完步驟一和步驟二后可以發(fā)現(xiàn),在執(zhí)行多脈沖平均的過程中,可能會(huì)存在多個(gè)活動(dòng)的平均窗口,如果下一個(gè)讀取的脈沖高度正好在包含平均和的兩個(gè)通道之間,則會(huì)隨機(jī)將脈沖高度添加到兩個(gè)通道之一; 如果下一個(gè)讀取的脈沖幅度所產(chǎn)生的平均窗口覆蓋了該區(qū)域內(nèi)已經(jīng)存在的某個(gè)脈沖幅度和,那么則用新的脈沖幅度與原有脈沖幅度相加,并將該平均窗口內(nèi)的脈沖個(gè)數(shù)加一。
如圖4所示,在讀幅度為659的脈沖之前,所示脈沖幅度和的區(qū)域內(nèi)已經(jīng)有650和660兩個(gè)脈沖幅度和,每一個(gè)脈沖幅度和對(duì)應(yīng)的脈沖數(shù)量都為一,還未達(dá)到窗口內(nèi)取平均的條件。下一個(gè)讀取的脈沖幅度659產(chǎn)生的平均窗口范圍為657~661,剛好包含了原本存在660,因此更新后的脈沖幅度和就為1 319,并且該窗口內(nèi)脈沖個(gè)數(shù)為2,達(dá)到了取平均的條件。通過脈沖幅度和除以脈沖個(gè)數(shù)計(jì)算出平均脈沖幅度為660,已經(jīng)被平均的脈沖幅度和以及脈沖個(gè)數(shù)清零,其余窗口內(nèi)的計(jì)數(shù)保持不變,同時(shí)將取平均值得到的脈沖幅度所對(duì)應(yīng)的道址上的計(jì)數(shù)加一,如圖5所示。
圖3 MPLA算法實(shí)現(xiàn)示例步驟二Fig.3 Step 2 of MPLA algorithm implementation example
圖4 MPLA算法實(shí)現(xiàn)示例步驟三Fig.4 Step 3 of MPLA algorithm implementation example
圖5 MPLA算法實(shí)現(xiàn)示例步驟四Fig.5 Step 4 of MPLA algorithm implementation example
理論上來說,傳統(tǒng)的多道譜(MCA)體現(xiàn)的是能量上的概率密度(pdf),每一個(gè)脈沖幅度都對(duì)最終得到的譜圖的計(jì)數(shù)率有所貢獻(xiàn),而MPLA將這種概率密度作為輸入,并通過在指定范圍內(nèi)的平均將其轉(zhuǎn)換為相同能量上新的概率密度。
設(shè)r為窗口平均參數(shù),n為平均參數(shù)的脈沖高度數(shù)。X0是一個(gè)通過多脈沖幅度和本地取平均值得到脈沖幅度,且X0的范圍為[X0-r,X0+r]。令概率密度函數(shù)為f(x),累積分布函數(shù)為F(x)。若X1是一個(gè)隨機(jī)變量描述的落在數(shù)理統(tǒng)計(jì)范圍內(nèi)的第二脈沖幅度,則X1與X0有相同的概率密度函數(shù),但區(qū)間為[X1-r,X1+r]。道址上的計(jì)數(shù)將在采集到2個(gè)數(shù)理統(tǒng)計(jì)范圍內(nèi)的脈沖幅度時(shí)進(jìn)行遞增,相應(yīng)的隨機(jī)變量P0=(X0+X1)/2。輸入概率密度函數(shù)f(x)將完全指定X0和X1的分布,以及MPLA轉(zhuǎn)換產(chǎn)生的P0的概率密度函數(shù)。
因?yàn)閒(x)是關(guān)于μ對(duì)稱的,所以對(duì)于任意值b可以得出f(μ-b)=f(μ+b)。同理,如果f(x)是關(guān)于μ對(duì)稱的,那么fp(y)也關(guān)于μ對(duì)稱。如果f(x)是對(duì)稱的,并且不斷的趨于平均值μ,那么P0的方差小于X0。推導(dǎo)過程如下:
在不損失一般性的前提下,我們假設(shè)為μ=0。該理論簡(jiǎn)化后為式(1)所示。
(1)
表達(dá)式(1)的左側(cè)可進(jìn)一步展開為式(2)所示。
(2)
(3)
由于f(x)是對(duì)稱的,并且一直朝著均值遞增(此處我們假定了均值為0),我們可以得出假如x0=μ=0那么E[X1|X0=x0]=x0,同樣的,假如x0≠μ那么|E[X1|X0=x0]|<|x0|,在此基礎(chǔ)上,式(3)可進(jìn)一步推導(dǎo)出式(4)。
(4)
上述推論表明,不管MPLA的參數(shù)值如何,MPLA轉(zhuǎn)換都不會(huì)改變分布,并保持對(duì)稱性。如果原始分布是單個(gè)高斯峰,則所得分布將具有相同的平均值,維持峰值位置,并且也將是對(duì)稱的。理想的MPLA轉(zhuǎn)換降低了原始分布的半高寬,銳化了峰形。對(duì)于單個(gè)高斯峰,半高寬與分布的方差有關(guān)。
MPLA算法是一種在線實(shí)時(shí)進(jìn)行的脈沖處理方法,其實(shí)施平臺(tái)在硬件上包括前端的探測(cè)器、后端的PC處理軟件以及核心部分的數(shù)字信號(hào)處理器。如前文所述,MPLA算法本身并沒有使用范圍的限制,它是將傳統(tǒng)的概率密度作為輸入,并通過在指定范圍內(nèi)的平均將其轉(zhuǎn)換為相同能量上新的概率密度,其本質(zhì)上就是一種概率密度的轉(zhuǎn)換,并沒有對(duì)探測(cè)器有任何的限定。本文所采用的探測(cè)器是AMPTEK生產(chǎn)的高性能FAST-SDD探測(cè)器。數(shù)字信號(hào)處理器包括信號(hào)轉(zhuǎn)換、多級(jí)放大、模數(shù)轉(zhuǎn)換等前端處理單元、脈沖幅度分析單元、MPLA單元以及脈沖成譜單元,其結(jié)構(gòu)框圖如圖6所示。
圖6 數(shù)字信號(hào)處理器結(jié)構(gòu)圖Fig.6 Structure chart of the digital signal processor
數(shù)字信號(hào)處理器的前端處理部分不是本文講述的重點(diǎn),因此下文中將重點(diǎn)介紹MPLA單元的硬件實(shí)現(xiàn)過程。MPLA單元在FPGA中實(shí)現(xiàn),采用的芯片是Xilinx公司Spartan-3家族的XC3S400-4TQG144C。如圖6所示,MPLA單元是整個(gè)數(shù)字信號(hào)處理器的核心部分,其算法的實(shí)現(xiàn)主要包含了算法控制、平均窗口定位、窗口內(nèi)平均以及相應(yīng)道址上的計(jì)數(shù)更新四個(gè)部分。具體的實(shí)現(xiàn)流程在前文中已經(jīng)有詳細(xì)的描述,此處不再贅述。在MPLA單元處理完成后,更新的計(jì)數(shù)則通過SPI總線傳遞給MCU,MCU與上位機(jī)之間通過CAN總線進(jìn)行通信。此外,MCU的另一SPI總線連接到前端處理電路中,用于控制前端電路中的放大參數(shù)和偏移量的調(diào)整。
由圖7可以看出經(jīng)過MPLA算法處理后的譜線半高寬明顯降低,也就是說能量分辨率有了明顯的提高,但與此同時(shí),多脈沖窄窗口求平均的成譜方法也造成了計(jì)數(shù)率的降低。
由于放射性衰變是隨機(jī)的,所以計(jì)數(shù)率總是在圍繞一個(gè)平均值不停地?cái)[動(dòng),即放射源在每單位時(shí)間內(nèi)發(fā)生衰變的原子數(shù)目是不同的,時(shí)多時(shí)少,有起有伏,但是它比較集中地在某一范圍內(nèi)波動(dòng),而這種現(xiàn)象就是放射性衰變的統(tǒng)計(jì)漲落。由概率統(tǒng)計(jì)理論可知,當(dāng)統(tǒng)計(jì)平均值較大時(shí),核衰變是服從正態(tài)分布的,這也是核衰變的重要統(tǒng)計(jì)特點(diǎn)。在實(shí)驗(yàn)環(huán)節(jié)中,以鐵礦樣品為測(cè)量對(duì)象,采用MPLA算法并進(jìn)行快速多脈沖成譜處理后的結(jié)果如圖7所示。
圖7 采用MPLA算法前后得到的譜圖Fig.7 Spectrum obtained before and after MPLA
圖中黑色譜線是通過傳統(tǒng)成譜方法得到的,紅色譜線是采用MPLA算法處理后得到的,相比之下,經(jīng)MPLA算法處理后的譜線有效地降低了特征峰的半高寬。這就驗(yàn)證了前文理論推導(dǎo)部分得出的結(jié)論,理想的MPLA轉(zhuǎn)換降低了原始分布的半高寬,銳化了峰形。
給出了多脈沖局部平均變換的理論描述及實(shí)現(xiàn)原理,證明了MPLA算法的一些基本特性,對(duì)具有單個(gè)峰的對(duì)稱分布,MPLA變換總是保持均值不變并減小FWHM。理論推導(dǎo)和實(shí)驗(yàn)結(jié)果表明,MPLA算法具有降低譜峰的FWHM的優(yōu)點(diǎn)。但與此同時(shí),MPLA的一個(gè)缺點(diǎn)是與傳統(tǒng)成譜方法相比,多脈沖平均的過程中使得計(jì)數(shù)減少。MPLA譜中計(jì)數(shù)減少的總數(shù)受MPLA算法中平均窗口內(nèi)脈沖幅度的個(gè)數(shù)控制。有很多脈沖處理技術(shù)可以彌補(bǔ)這種計(jì)數(shù)損失,例如,可以收集大量的脈沖高度測(cè)量,然后重復(fù)、隨機(jī)地進(jìn)行采樣以應(yīng)用MPLA變換。筆者將在接下來的工作中對(duì)這些脈沖處理技術(shù)進(jìn)行更深入的研究,以期為MPLA算法找到一種非常契合的脈沖處理技術(shù)來補(bǔ)償計(jì)數(shù)率損失。即便存在計(jì)數(shù)損失的缺陷,但不可否認(rèn),MPLA算法為降低探測(cè)器噪聲提供了一種實(shí)時(shí)、高效、通用的處理方法,并且具有可證明的理論保證,這對(duì)提高探測(cè)器的能量分辨率具有極大的應(yīng)用意義。