本發(fā)明涉及醫(yī)療器械的信號處理領域，更具體地涉及數(shù)字pet領域的一種閃爍脈沖數(shù)字化信號的擬合方法。
背景技術：
正電子發(fā)射斷層成像(positronemissiontomography，全文均簡稱pet)通過捕捉人體內(nèi)因正電子湮滅而發(fā)出的γ光子獲取以正電子核素為標記的示蹤劑在人體內(nèi)的分布情況，進而獲取臟器功能，代謝等病理生理特征。獲取γ光子所攜帶的能量、位置以及時間信息的準確性直接影響到系統(tǒng)成像的性能。閃爍探測器用于高能粒子(例如γ、x等粒子)的捕獲及相應的能量沉積信息測量，由于其具有探測效率高、時間響應快和信息測量準確等特點，廣泛應用于輻射探測與成像領域。閃爍探測器通常由閃爍體(例如bgo、lyso、labr等)和光電轉(zhuǎn)換器(例如pmt、sipm、apd等)構(gòu)成，入射的γ、x等粒子與閃爍體相互作用產(chǎn)生熒光，光電轉(zhuǎn)換器件則將熒光轉(zhuǎn)換為相應的閃爍脈沖。將閃爍脈沖數(shù)字化，利用數(shù)字信號處理技術分析數(shù)字閃爍脈沖以提取粒子的能量沉積信息，通過進一步的分析數(shù)字閃爍脈沖、優(yōu)化和升級處理算法即可不斷提高粒子能量沉積信息提取的精度，提高閃爍探測器的性能，進而提高pet系統(tǒng)的成像質(zhì)量。
數(shù)字閃爍脈沖采樣的一種方法是多閾值采樣(multi-voltagethreshold，以下簡稱mvt)方法，該方法根據(jù)閃爍脈沖的特征合理的設置多個電壓閾值，通過對閃爍脈沖越過電壓閾值的時間的數(shù)字化最終實現(xiàn)閃爍脈沖的數(shù)字化采樣。利用少量的閾值比較器配合時間數(shù)字轉(zhuǎn)換器(time-to-digitalconverter，簡稱tdc)便可設計出相應的mvt采樣電路進行閃爍脈沖的mvt采樣，完成相應的數(shù)字閃爍脈沖的獲取。mvt方法用于閃爍脈沖的數(shù)字化和粒子能量沉積信息的提取，具有性能優(yōu)異、工程實現(xiàn)簡單、硬件成本低廉等優(yōu)點。現(xiàn)有研究結(jié)果表明，采用3～8個閾值的mvt采樣電路獲取的數(shù)字閃爍脈沖，配合基于先驗的閃爍脈沖形狀信息的數(shù)字脈沖分析方法即可實現(xiàn)粒子能量沉積信息的精確獲取。
現(xiàn)有的數(shù)字脈沖分析方法之一為最小二乘算法，其包括牛頓法、梯度法等。然而，現(xiàn)有的大部分最小二乘算法在對pet系統(tǒng)獲取到的數(shù)字閃爍脈沖信號進行擬合時耗時極長，加上大數(shù)量級的閃爍脈沖數(shù)目，導致實際處理過程中效率低下，從采集數(shù)據(jù)到得到醫(yī)學圖像的耗時過長，已經(jīng)逐漸無法滿足工業(yè)需求，而且這些算法具有一定的使用限定，無法在穩(wěn)定性、精確性上達到令人滿意的平衡。另外有一些擬合算法是運用最小二乘的思想，通過選定雙指數(shù)模型中的兩個參數(shù)來進行擬合，以犧牲脈沖完整性的代價來實現(xiàn)一定的時間效率，但是，所獲取的閃爍脈沖時間、能量信息并不十分精確，最終在pet系統(tǒng)的時間分辨率和能量分辨率上會有一定惡化。
列文伯格-馬夸爾特(levenberg-marquardt)法是最優(yōu)化算法中的一種，也是使用最廣泛的非線性最小二乘迭代算法，它是利用梯度求最大(小)值，介于牛頓法與梯度下降法之間的一種非線性優(yōu)化方法，并且同時具有梯度法和牛頓法的優(yōu)點，因此，它的應用領域非常廣泛，如：經(jīng)濟學、管理優(yōu)化、網(wǎng)絡分析、最優(yōu)設計、機械或電子設計等等。因此，為了進一步提高脈沖重建后的圖像質(zhì)量，如果能夠?qū)⒘形牟?馬夸爾特算法應用于經(jīng)mvt采樣獲得的閃爍脈沖數(shù)字化信號的處理，將是兼顧效率、準確度、穩(wěn)定性等綜合因素的最佳擬合方法，能夠更好的滿足工業(yè)需求。
技術實現(xiàn)要素：
本發(fā)明的目的是提供一種閃爍脈沖數(shù)字化信號的擬合方法，從而解決現(xiàn)有技術的數(shù)字脈沖擬合方法不能兼顧效率、準確度和穩(wěn)定性的問題。
為了解決上述技術問題，本發(fā)明的技術方案是提供一種閃爍脈沖數(shù)字化信號的擬合方法，該擬合方法包括以下步驟：
步驟s1：選定閃爍脈沖的先驗模型為雙指數(shù)模型v(t)＝a·ebt+c·edt，將所述雙指數(shù)模型確定為目標函數(shù)；
步驟s2：選定閃爍脈沖數(shù)據(jù)庫中的經(jīng)典脈沖的采樣點，所述采樣點包括八個電壓閾值序列和八個時間點序列，其中，八個所述電壓閾值序列為v1、v2、v3、v4、v5、v6、v7和v8，八個所述時間點序列為t1、t2、t3、t4、t5、t6、t7和t8；
步驟s3：將所述電壓閾值序列存為電壓數(shù)組v＝[v1，v2，v3，v4，v5，v6，v7，v8]，對所述雙指數(shù)模型的曲線進行時間上的平移，所述平移的具體步驟為：將第一個時間點t1置0，后續(xù)時間點分別與所述第一個時間點t1求差，然后將得到的數(shù)據(jù)存為時間數(shù)組t＝[0，(t2-t1)，(t3-t1)，(t4-t1)，(t5-t1)，(t6-t1)，(t7-t1)，(t8-t1)]；
步驟s4：將所述目標函數(shù)的四個初始默認參數(shù)a、b、c、d設置為1，存為初始參數(shù)數(shù)組para＝[1，1，1，1]；
步驟s5：將所述電壓數(shù)組v、時間數(shù)組t、所述脈沖采樣點個數(shù)n、所述初始參數(shù)數(shù)組para、所述目標函數(shù)v(t)輸入到列文伯格-馬夸爾特擬合函數(shù)中，進行擬合得到四個擬合參數(shù)a’、b’、c’、d’，將四個所述擬合參數(shù)存為第二參數(shù)數(shù)組paranew＝[a’，b’，c’，d’]；
步驟s6：對于閃爍脈沖數(shù)據(jù)庫中的每個閃爍脈沖，除了步驟s4中用所述第二參數(shù)數(shù)組paranew代替所述初始參數(shù)數(shù)組para外，依次重復所述步驟s1-步驟s5，直至完成所有閃爍脈沖的擬合，記第k個閃爍脈沖擬合后的函數(shù)為
步驟s7：對擬合后確定擬合參數(shù)的閃爍脈沖的目標函數(shù)vk(t)積分得到第k個脈沖的能量值ek；
步驟s8：解方程vk(t)＝0得到曲線過零點tk0，得到第k個脈沖的到達時間tk＝t1+tk0，脈沖的時間信息、能量信息提取完畢。
閃爍脈沖的先驗模型通過數(shù)字示波器采集得到，閃爍脈沖的曲線包括快速的上升沿和緩慢的下降沿。
閃爍脈沖數(shù)據(jù)庫通過閃爍脈沖獲取平臺記錄數(shù)字閃爍脈沖形成。
閃爍脈沖獲取平臺包括一對相對設置的閃爍探測器以及一個射源，射源置于一對所述閃爍探測器軸向連線的中間位置，所述閃爍探測器分別通過同軸電纜接入所述數(shù)字示波器。
閃爍探測器為si-bdm探測器，所述射源為點狀18fdg射源。
數(shù)字示波器以16ghz的模擬帶寬和50gsps的采樣率進行數(shù)字化采樣，記錄采集到的閃爍脈沖。
列文伯格-馬夸爾特擬合函數(shù)的原型為voidlmcurve(intn_par，double*par，intm_dat，constdouble*t，constdouble*y，double(*f)(doublet，constdouble*par)，constlm_control_struct*control，lm_status_struct*status)，其中，intn_par表示待擬合模型函數(shù)的參數(shù)個數(shù)；double*par表示參數(shù)初始值；intm_dat表示擬合采樣點的個數(shù)；constdouble*t表示采樣點的橫坐標數(shù)組；constdouble*y表示采樣點的縱坐標數(shù)組；double(*f)(doublet，constdouble*par)表示擬合依據(jù)的曲線模型；constlm_control_struct*control表示擬合算法的控制參數(shù)；lm_status_struct*status表示擬合算法的狀態(tài)參數(shù)。
數(shù)據(jù)類型為包括double和float，對于性能好、配置高的服務器，數(shù)據(jù)類型設置為double，否則數(shù)據(jù)類型設置為float。
目標函數(shù)vk(t)的積分范圍為0-200。
曲線過零點tk0＝[1/(dk-bk)*ln(-ak/ck)]。
本發(fā)明的閃爍脈沖數(shù)字化信號的擬合方法，利用列文伯格-馬夸爾特法對基于mvt方法獲取的數(shù)字閃爍脈沖信號進行非線性擬合，選用雙指數(shù)模型，通過從已有的脈沖數(shù)據(jù)庫中提取的雙指數(shù)模型參數(shù)作為算法的初值，能大大提高算法運行時間，同時算法本身易于openmp/gpu進行加速，隨著計算機的更新?lián)Q代算法本身的時間性能也會逐漸提高。通過本方法對脈沖進行符合并重建后，得到的能量分辨率有顯著提升，本方法不僅能夠提高處理效率，而且通過本方法還能保證信息提取的準確度和穩(wěn)定性，因此，本發(fā)明是在效率、準確、穩(wěn)定等綜合因素下的最佳擬合方法。
附圖說明
圖1是根據(jù)本發(fā)明的閃爍脈沖數(shù)字化信號的擬合方法的脈沖的先驗模型示意圖；
圖2是根據(jù)本發(fā)明一個實施例的閃爍脈沖獲取平臺的布置示意圖；
圖3是根據(jù)本發(fā)明一個實施例的經(jīng)典閃爍脈沖的示意圖；
圖4是根據(jù)本發(fā)明一個實施例的采用雙指數(shù)模型的閃爍脈沖數(shù)字化信號的擬合方法應用于pet信號處理的效果圖；
圖5是現(xiàn)有技術的采用直線指數(shù)模型所獲得的能量分辨率的效果圖。
具體實施方式
以下結(jié)合具體實施例，對本發(fā)明做進一步說明。應理解，以下實施例僅用于說明本發(fā)明而非用于限制本發(fā)明的范圍。
本發(fā)明提供的閃爍脈沖數(shù)字化信號的擬合方法，包括以下步驟：
步驟s1：選定閃爍脈沖的先驗模型為雙指數(shù)模型v(t)＝a·ebt+c·edt，將該雙指數(shù)模型確定為目標函數(shù)；
步驟s2：選定閃爍脈沖數(shù)據(jù)庫中的經(jīng)典脈沖的采樣點，該采樣點包括八個電壓閾值序列和八個時間點序列，其中，八個電壓閾值序列為v1、v2、v3、v4、v5、v6、v7和v8，八個時間點序列為t1、t2、t3、t4、t5、t6、t7和t8；
步驟s3：將上述電壓閾值序列存為電壓數(shù)組v＝[v1，v2，v3，v4，v5，v6，v7，v8]，對雙指數(shù)模型的曲線進行時間上的平移，平移的具體方法為：將第一個時間點t1置0，后續(xù)時間點分別與t1求差，將得到的數(shù)據(jù)存為時間數(shù)組t＝[0，(t2-t1)，(t3-t1)，(t4-t1)，(t5-t1)，(t6-t1)，(t7-t1)，(t8-t1)]；
步驟s4：選定目標函數(shù)的四個初始默認參數(shù)a、b、c、d為1，存為初始參數(shù)數(shù)組para＝[1，1，1，1]；
步驟s5：將電壓數(shù)組v、時間數(shù)組t、脈沖采樣點個數(shù)n(n＝8)、初始參數(shù)數(shù)組para、目標函數(shù)v(t)輸入到擬合函數(shù)lmcurve中，根據(jù)需要設置相應的數(shù)據(jù)類型、擬合速度等參數(shù)，進行擬合，得到的輸出結(jié)果為四個新的擬合參數(shù)a’、b’、c’、d’，將這四個新的擬合參數(shù)存為第二參數(shù)數(shù)組paranew＝[a’，b’，c’，d’]；
步驟s6：針對實際應用中得到的多個閃爍脈沖數(shù)據(jù)的擬合過程，除了步驟s4中用第二參數(shù)數(shù)組paranew來代替初始參數(shù)數(shù)組para外，其余的閃爍脈沖處理依次重復步驟s1-步驟s5，直至完成所有閃爍脈沖的擬合，記第k個閃爍脈沖擬合后的函數(shù)為
步驟s7：對擬合后確定擬合參數(shù)的閃爍脈沖的目標函數(shù)vk(t)積分得到第k個脈沖的能量值ek；
步驟s8：解方程vk(t)＝0得到曲線過零點tk0(為負數(shù))，得到第k個脈沖的到達時間tk＝t1+tk0，脈沖的時間信息、能量信息提取完畢。
上述步驟s1中，脈沖的先驗模型通過數(shù)字示波器采集得到，如圖2所示，該先驗模型的閃爍脈沖曲線包括快速的上升沿和緩慢的下降沿，該閃爍脈沖的形狀信息可通過直線指數(shù)和雙指數(shù)模型描述。本發(fā)明中采用了雙指數(shù)模型進行描述，由于閃爍脈沖的下降沿非常快，雙指數(shù)模型能有效地提高時間分辨率計算的準確性。
上述步驟s2中，閃爍脈沖數(shù)據(jù)庫通過閃爍脈沖獲取平臺記錄數(shù)字閃爍脈沖形成。圖2為根據(jù)本發(fā)明一個實施例的閃爍脈沖獲取平臺1的布置示意圖，該閃爍脈沖獲取平臺1包括一對相對設置的閃爍探測器10以及一個射源20，在圖2的實施例中，所采用的閃爍探測器10為si-bdm探測器，射源20為點狀18fdg射源，該點狀18fdg射源20置于一對閃爍探測器10軸向連線的中間位置；兩個閃爍探測器10分別通過同軸電纜40接入數(shù)字示波器30，比如，同軸電纜40分別接入數(shù)字示波器的通道二和通道三；數(shù)字示波器30以16ghz的模擬帶寬和50gsps的采樣率進行數(shù)字化采樣，記錄通過兩個通道采集到的閃爍脈沖，從而形成閃爍脈沖數(shù)據(jù)庫。
根據(jù)本發(fā)明的一個實施例，上述步驟s2中采用的經(jīng)典脈沖及采樣點的設置如圖3所示，其中，八個電壓閾值序列為v1、v2、v3、v4、v5、v6、v7和v8，八個時間點序列為t1、t2、t3、t4、t5、t6、t7和t8，圖3中電壓閾值的具體值的設定可以根據(jù)實際需要進行調(diào)整，在此不在贅述。
上述步驟s3中，數(shù)據(jù)的存儲可通過適當?shù)牟僮髌脚_和編程語言實現(xiàn)，比如c++語言，在此不再贅述。
上述步驟s5中，選用的列文伯格-馬夸爾特(levenberg-marquardt)擬合算法的函數(shù)原型為voidlmcurve(intn_par，double*par，intm_dat，constdouble*t，constdouble*y，double(*f)(doublet，constdouble*par)，constlm_control_struct*control，lm_status_struct*status)，該函數(shù)原型的參數(shù)中，intn_par表示待擬合模型函數(shù)的參數(shù)個數(shù)；double*par表示參數(shù)初始值；intm_dat表示擬合采樣點的個數(shù)；constdouble*t表示采樣點的橫坐標數(shù)組；constdouble*y表示采樣點的縱坐標數(shù)組；double(*f)(doublet，constdouble*par)表示擬合依據(jù)的曲線模型；constlm_control_struct*control表示擬合算法的控制參數(shù)，比如擬合速度、容錯等參數(shù)，一般取默認值；lm_status_struct*status表示擬合算法的狀態(tài)參數(shù)。
上述步驟s5中，根據(jù)需要設置相應的數(shù)據(jù)類型、擬合速度等參數(shù)是指根據(jù)所需擬合的閃爍脈沖的數(shù)目大小以及系統(tǒng)服務器的內(nèi)存配置，對于性能好、配置高的服務器，可設置數(shù)據(jù)類型為double，否則設置數(shù)據(jù)類型為float；同時，該擬合算法自帶的控制參數(shù)lm_control_struct*control中包含有擬合速度、步階、擬合相對誤差等參數(shù)，這些參數(shù)都可以根據(jù)當前系統(tǒng)服務器、硬件配置等進行調(diào)試，以實現(xiàn)提高效率和準確性的目的，例如，對于本發(fā)明的一個實施例，可設置數(shù)據(jù)類型為double，其余的擬合算法的控制參數(shù)設置為默認值。
上述步驟s7中，對vk(t)積分的范圍0-200。
上述步驟s8中，解方程是令
圖4為根據(jù)本方法對一個閃爍脈沖進行數(shù)據(jù)處理后得到的能量分辨率圖譜，圖5為利用直線指數(shù)模型獲得的閃爍脈沖的能量分辨率的效果圖，對比圖4和圖5可以看出，本發(fā)明采用雙指數(shù)模型得到的閃爍脈沖的能量分辨率為17.9％，相比于采用直線指數(shù)模型的能量分辨率23.4％，顯然本發(fā)明獲得的能量分辨率具有顯著提升，這對數(shù)字pet系統(tǒng)而言是明顯的進步。
本發(fā)明基于mvt方法得到閃爍脈沖數(shù)據(jù)庫，根據(jù)預設的四個電壓閾值，在脈沖上升沿和下降沿都得到四個采樣點，即一個脈沖用八個采樣點(時間，電壓序列)表示，通過將已有的脈沖數(shù)據(jù)庫中提取的雙指數(shù)模型參數(shù)作為擬合算法的初值，能大大提高算法運行時間，該方法本身也易于openmp/gpu進行加速，該方法做到了穩(wěn)定性和精確性的平衡，在完整還原雙指數(shù)模型的閃爍脈沖情況下提高了效率和算法穩(wěn)定性，根據(jù)此方法最后得到的pet系統(tǒng)能量分辨率有顯著提高。
以上所述的，僅為本發(fā)明的較佳實施例，并非用以限定本發(fā)明的范圍，本發(fā)明的上述實施例還可以做出各種變化。即凡是依據(jù)本發(fā)明申請的權(quán)利要求書及說明書內(nèi)容所作的簡單、等效變化與修飾，皆落入本發(fā)明專利的權(quán)利要求保護范圍。本發(fā)明未詳盡描述的均為常規(guī)技術內(nèi)容。