專利名稱：：奇異譜分析截?cái)囝l譜信號(hào)重建方法及其應(yīng)用的制作方法
技術(shù)領(lǐng)域：
：本發(fā)明屬于截?cái)囝l譜數(shù)字信號(hào)重建技術(shù)，特別用于磁共振截?cái)囝l譜成象。儀器儀表在采集或傳輸信號(hào)時(shí)，由于儀器儀表或傳輸系統(tǒng)的頻寬有限，或采集條件苛刻等，使得采集到的信號(hào)缺乏高頻成分，只能獲取低頻分量的信號(hào)?；蛘攉@得足夠高的高頻成分要化費(fèi)很大的設(shè)備開消(儀器儀表的生產(chǎn)成本或傳輸設(shè)備成本很高)或時(shí)間開消(信號(hào)采集時(shí)間很長(zhǎng))。這就出現(xiàn)了截?cái)囝l譜信號(hào)(缺乏高頻成分的信號(hào))的重建問題，或者說從已有的截?cái)囝l譜數(shù)據(jù)中重建出完整無缺頻譜數(shù)據(jù)的問題。例如磁共振成象(MRI)技術(shù)是一種日益令人注目的醫(yī)學(xué)圖象診斷手段。但其因成象時(shí)間長(zhǎng)，信噪比很難提高，而導(dǎo)致MRI設(shè)備十分昂貴，嚴(yán)重制約著它的廣泛應(yīng)用。為了縮小采集信號(hào)時(shí)間、提高信噪比。人們?cè)噲D用磁共振截?cái)囝l譜代替完全頻譜成象。其方法是只采集低頻頻譜數(shù)據(jù)，其余高頻頻譜數(shù)據(jù)用零填補(bǔ)，然后用付里葉反變換成象(例如深圳安科公司生產(chǎn)的MRI設(shè)備)。這種截?cái)囝l譜成象方法叫補(bǔ)零法，其缺陷是圖象中有偽影、圖象質(zhì)量差。這偽影由于截?cái)囝l譜欠缺高頻頻譜分量所造成的，所以國(guó)際醫(yī)學(xué)成象界把這種偽影叫截?cái)鄠斡?。由磁共振截?cái)囝l譜重建出與完全頻譜重建一樣的磁共振圖象，就可以既不增加硬件資源，又不損失任何圖象信息的情況下，節(jié)約掃描時(shí)間，提高磁共振設(shè)備輸出信號(hào)的信噪比，從而降低設(shè)備性能要求，降低設(shè)備成本。目前國(guó)內(nèi)外尚無一種既能消除偽影又能滿足臨床使用要求的成象方法。解決磁共振截?cái)囝l譜成象中的截?cái)鄠斡皢栴}是磁共振截?cái)囝l譜成象中急需解決的關(guān)鍵問題，是目前國(guó)際磁共振醫(yī)學(xué)成象界的主要研究熱點(diǎn)之一，是一項(xiàng)有重要學(xué)術(shù)價(jià)值、能產(chǎn)生巨大經(jīng)濟(jì)效益和促進(jìn)社會(huì)醫(yī)療保健事業(yè)發(fā)展的高難度新技術(shù)。本發(fā)明的目的是針對(duì)現(xiàn)有技術(shù)的問題，從截?cái)囝l譜中重建出完整頻譜，從而解決截?cái)囝l譜信號(hào)重建中的偽影問題；應(yīng)用本方法于磁共振截?cái)囝l譜成象領(lǐng)域，效果顯著并能降低設(shè)備成本。為實(shí)現(xiàn)上述目的，采用以下發(fā)明構(gòu)思技術(shù)關(guān)鍵是按數(shù)學(xué)理論中，函數(shù)可表示為某個(gè)函數(shù)系的函數(shù)的泛函的原理，建立信號(hào)譜奇異譜分析數(shù)學(xué)模型，從采集到的已有的截?cái)囝l譜數(shù)據(jù)中抽取模型參數(shù)(奇異點(diǎn)和奇異譜分析的加權(quán)系數(shù))，再根據(jù)奇異譜分析數(shù)學(xué)模型，推出截?cái)囝l譜所隱含的高頻分量，解決截?cái)囝l譜信號(hào)重建中的偽影問題。發(fā)明的主要內(nèi)容包括1、奇異譜分析數(shù)學(xué)模型我們知道，在一定條件下，一個(gè)函數(shù)可表示為三角函數(shù)系、小波函數(shù)系或冪函數(shù)系的線性泛函。即在一定條件下一個(gè)函數(shù)可用另一函數(shù)序列的加權(quán)和表示。對(duì)于任一離散頻譜函數(shù)，能否用一個(gè)離散函數(shù)序列的線性泛函表示？回答是肯定的。為此，我們給出如下定義定義1、定義離散函數(shù)空間{Wk(i)，k，i為整數(shù)}為奇異空間，如果它滿足以下三個(gè)條件(1)Wk(i)是Wo(i)的平移或彭脹k生成的；(2)Wk(i)只有一個(gè)唯一的差分不為零的點(diǎn)k(奇異點(diǎn)k)(3)Wk(i)的Fourier變換Wk^(u)=F[Wk(i)]]]>存在。(F[.]為Fourier變換算符)并稱Wk(i)為奇異函數(shù)，為奇異譜函數(shù)，函數(shù)集合{Wk^(u)|k,u]]>為整數(shù)}為奇異譜空間。按上述定義，數(shù)學(xué)上容易得出下面結(jié)論1、奇異空間{Wk(i)，k，i為整數(shù)}是完備的，即任何離散信號(hào)f(j)，j＝0，1，...，N-1可表示空間{Wk(i)}中的線性泛函。即f(j)=Σk=0N-1akwk(j),j∈{0,1,...,N-1}......(1)]]>其中式(1)的各加權(quán)系數(shù)為ak={f(k)-f(k-1)k=1,2,...,N-1f(k)k=0;......(2)]]>2、任何一個(gè)離散函數(shù)f(j)的頻譜函數(shù)F(u)都可以表示為奇異譜空間{Wk^(u)|k,]]>u為整數(shù)}中的線性泛函。即F(u)=Σk=0N-1akWk^(u),u∈{0,1,...,N-1}......(3)]]>3、若信號(hào)中僅有Q個(gè)奇異點(diǎn)(在ak，k＝0、1，...，N-1中僅有Q個(gè)不為0)，則信號(hào)的譜函數(shù)F(u)可由{Wk^(u)|k=0,1,...,N-1}]]>中對(duì)應(yīng)的Q個(gè)奇異譜函數(shù)的線性泛函表示。即F(u)=Σi=1QabiWbi^(u),abi≠0......(4)]]>其中是以bi位置為奇異點(diǎn)的奇異譜函數(shù)，abi是的權(quán)系數(shù)。所以，任何信號(hào)的頻譜可以由信號(hào)的奇異點(diǎn)及其權(quán)系數(shù)唯一確定。2、奇異譜分析截?cái)囝l譜信號(hào)重建的算法按照信號(hào)奇異譜分析數(shù)學(xué)模型的第3個(gè)結(jié)論知，獲得信號(hào)的奇異點(diǎn)及其奇異譜分析加權(quán)系數(shù)就等價(jià)于獲取信號(hào)的頻譜數(shù)據(jù)。如果能從磁共振截?cái)囝l譜中析取信號(hào)的奇異點(diǎn)及奇異譜分析加權(quán)系數(shù)，那么就能由截?cái)囝l譜重建磁共振頻譜的高頻頻譜分量，從而解決磁共振成象中存在偽影問題。我們的思路是先選用恰當(dāng)?shù)恼恍〔ê瘮?shù)，對(duì)已有的截?cái)囝l譜求小波變換系數(shù)；再通過調(diào)整尺度參數(shù)，測(cè)定小波系數(shù)(變換值)的極點(diǎn)，即信號(hào)的奇異點(diǎn)；最后由已知磁共振截?cái)囝l譜，按結(jié)論3確定信號(hào)的奇異譜分析加權(quán)系數(shù)。最后用FOURIER反變換重建磁共振圖象。(1)奇異點(diǎn)測(cè)定對(duì)于N維離散信號(hào)情況，不難證明對(duì)于頻譜函數(shù)為的信號(hào)的小波系數(shù)cj，k可以寫成cj,k=1NΣu=-N/2N/2f^(u)ψ^j,k(u)‾---(5)]]>我們可以通過信號(hào)的頻譜求出信號(hào)的小波系數(shù)。為正交二進(jìn)B小波(簡(jiǎn)稱B小波)的付里葉變換。由于信號(hào)能量大部分集結(jié)在低頻部分，所以由一定數(shù)量頻譜數(shù)據(jù)的低頻分量求出的小波系數(shù)可以用來確定象信號(hào)的奇異點(diǎn)，比如，用小于二分之一截止頻率的低頻分量。當(dāng)信號(hào)在某點(diǎn)岀現(xiàn)奇異點(diǎn)時(shí)，則該點(diǎn)將是小波變換系數(shù)的模的極值，且這模的極值將隨著尺度的增加而增加。若是噪聲引起的小波變換系數(shù)的模極值，將隨著尺度的增加而減小。據(jù)此就可確定奇異點(diǎn)，且受噪聲影響少。(2)權(quán)系數(shù)確定設(shè)N維信號(hào)f的截?cái)囝l譜為F=[f^(0),f^(1),...,f^(n-1)]T,]]>其中n＜N。由信號(hào)的截?cái)囝l譜測(cè)得的奇異點(diǎn)為bi，i＝1，2，...，Q。由奇異譜分析數(shù)學(xué)模型的結(jié)論3知，信號(hào)f的Fourier系數(shù)為f^(u)=Σi=1QybiWbi^(u),u∈{0,1,...,N-1}......(6)]]>其中ybi為待定權(quán)系數(shù)。顯然，對(duì)于截?cái)囝l譜F=[f^(0),f^(1),...,f^(n-1)]T]]>有方程組記y＝{yb1，yb2，...，ybQ}T為Q維權(quán)矢量，式(7)的矩陣形式為W^·y=F......(8)]]>其中是由奇異點(diǎn)確定的、已知的矩陣。顯然，在奇異點(diǎn)測(cè)定正確的前提下，方程組(8)是相容，并且當(dāng)n＞Q時(shí)，可以證明的秩為Q，按線性方程組解的性質(zhì)可得唯一解，高頻頻譜得到準(zhǔn)確恢復(fù)。當(dāng)n＜Q時(shí)，可以證明的秩為n，按線性方程組解的性質(zhì)可得唯一偽逆解，高頻頻譜得到最佳恢復(fù)。根據(jù)上述發(fā)明構(gòu)思設(shè)N維數(shù)字信號(hào)為的截?cái)囝l譜為F=[f^(0),f^(1),...,f^(n-1)]T]]>其中n＜N。則奇異譜分析截?cái)囝l譜信號(hào)重建方法的操作步驟為1、運(yùn)用下式計(jì)算各尺度下的小波系數(shù)cj,k=1NΣu=0n-1f^(u)ψ^j,k(u)‾,k=0,1...,N-1,]]>其中尺度j取0，1，2，32、檢索極值將隨著尺度的增加而增加的極點(diǎn)，即確定奇異點(diǎn)bi，i＝1，2，...，Q。3、用偽逆矩陣法求解以下方程組，求加權(quán)系數(shù)y＝{yb1，yb2，...，ybQ}T4、用下式計(jì)算截?cái)囝l譜的各高頻分量f^(u)Σi=1QybiWbi^(u),u=n,n+1,...,N-1;]]>5、對(duì)頻譜矢量F=[f^(0),f^(1),...,f^(n-1),f^(n),...,f^(N-1)]T]]>進(jìn)行Fourier反變換，得重建信號(hào)。設(shè)磁共振截?cái)囝l譜對(duì)應(yīng)的圖象尺寸為N×N，截?cái)囝l譜的尺寸為n×N，其中n＜Ne上述奇異譜分析截?cái)囝l譜信號(hào)重建方法應(yīng)用于磁共振截?cái)囝l譜成象的操作步驟為1、從第1行到第n行頻譜數(shù)據(jù)進(jìn)行Fourier反變換2、從第1列到第N列截?cái)囝l譜數(shù)據(jù)F=[f^(0),f^(1),...,f^(n-1)]T]]>反復(fù)執(zhí)行下例步驟(1)運(yùn)用上述方法重建截?cái)囝l譜F=[f^(0),f^(1),...,f^(n-1)]T]]>的信號(hào)x＝[x0，x1，...XN-1]T(2)求信號(hào)的模；3、輸出圖象。參照本發(fā)明的顯著進(jìn)步圖1表示無噪聲計(jì)算機(jī)模擬磁共振頻譜數(shù)據(jù)重建情況；圖2表示含5％高斯噪聲計(jì)算機(jī)模擬磁共振頻譜數(shù)據(jù)重建情況；圖3表示為低頻截?cái)囝l譜數(shù)據(jù)圖象。我們把奇異譜分析截?cái)囝l譜信號(hào)重建方法的進(jìn)行了大量實(shí)驗(yàn)。所用的兩類數(shù)據(jù)是實(shí)際磁共振截?cái)囝l譜數(shù)據(jù)和計(jì)算機(jī)模擬磁共振截?cái)囝l譜數(shù)據(jù)，計(jì)算機(jī)模擬的是人頭模型圖象(圖象尺寸是128×128，如圖1(a)所示)的截?cái)囝l譜數(shù)據(jù)。實(shí)際截?cái)囝l譜數(shù)據(jù)是由深圳安科公司提供的核磁共振人頭截面圖象(圖象尺寸是256×256)的截?cái)囝l譜數(shù)據(jù)。截?cái)囝l譜截取原頻譜數(shù)據(jù)的低頻部分，數(shù)據(jù)量是原總頻譜數(shù)據(jù)量的一半。為了充分認(rèn)識(shí)本發(fā)明的突出效果，我們用三種方式把它重建的圖象和補(bǔ)零法重建的圖象進(jìn)行了比較。第一是圖象直接比較第二是線圖比較(以圖象平面中某一直線上的灰度變化用一條曲線來表示，并用以這條直線上的位置作為橫坐標(biāo)，以這位置上的灰度作為縱坐標(biāo)，為方便我們稱之為線圖)第三是進(jìn)行保真計(jì)算。對(duì)于受某種污染的頻譜數(shù)據(jù)，以某種方法重建的圖象與原始標(biāo)準(zhǔn)圖象的接近程度，來表示這種圖象重建算法的優(yōu)劣。我們采用以下二種保真指標(biāo)，以區(qū)別圖象優(yōu)劣歸一化均方誤差，即σ=ΣI=1MΣJ=1N(xI,JO-xI,J)2ΣI=1MΣJ=1N(xI,JO-μ)2......(9)]]>歸一化絕對(duì)誤差，即E=ΣJ=1MΣJ=1N|xI,JO-xI,J|ΣI=1MΣJ=1N|xI,JO-μ|---(10)]]>其中XIJO、XI，J分別是標(biāo)準(zhǔn)圖象和重建圖象的第I行第J列的圖象素灰度值。μ是標(biāo)準(zhǔn)圖象的灰度均值。M是圖象的行數(shù)，N是圖象的列數(shù)。(1)計(jì)算機(jī)模擬截?cái)囝l譜數(shù)據(jù)成象的圖象比較無噪聲計(jì)算機(jī)模擬磁共振頻譜數(shù)據(jù)重建情況圖I(a)是完整頻譜重建的圖象，圖1(b)是截?cái)囝l譜用補(bǔ)零法重建的圖象，圖1(c)是本發(fā)明方法重建的圖象，圖象尺寸都是128×128。圖1(d)(e)和(f)分別是圖1(a)(b)和(c)圖的第42行的線圖。其截?cái)囝l譜數(shù)據(jù)為完整頻譜數(shù)據(jù)中的低頻64列。重建比較由圖1(a)(b)和(c)比較，圖1(d)(e)和(f)線圖比較得奇異譜成象法能完全消除截?cái)鄠斡埃杀?知奇異譜成象法的保真指標(biāo)比補(bǔ)零法高岀四個(gè)數(shù)量級(jí)。含5％高斯噪聲計(jì)算機(jī)模擬磁共振頻譜數(shù)據(jù)重建情況圖2(a)是完整頻譜重建的圖象，圖2(b)是截?cái)囝l譜用補(bǔ)零法重建的圖象，圖2(c)是本發(fā)明方法重建的圖象，圖象尺寸都是128×128。圖2(d)(e)和(f)分別是圖2(a)(b)和(c)圖的第42行的線圖。其截?cái)囝l譜數(shù)據(jù)為完整頻譜數(shù)據(jù)中的低頻64列。由圖2(a)(b)和(c)比較，圖2(d)(e)和(f)線圖比較得奇異譜成象法能完全消除截?cái)鄠斡?，由?知奇異譜成象法的保真指標(biāo)比補(bǔ)零法高出2個(gè)數(shù)量級(jí)。由表1、表2比較可見算法對(duì)噪聲敏感，這是因?yàn)橛捎谠肼暤囊胧沟闷娈慄c(diǎn)數(shù)量Q急劇上升，頻譜數(shù)據(jù)量n小于奇異點(diǎn)數(shù)量Q或方程(8)成為不相容，使方程(8)只能得一偽逆解。(2)實(shí)際磁共振截?cái)囝l譜數(shù)據(jù)成象的圖象比較圖3(a)(b)圖的截?cái)囝l譜數(shù)據(jù)為低頻80列，(a)為補(bǔ)零法重建，(b)為奇異譜分析法重建。截?cái)囝l譜數(shù)據(jù)直接來源于深圳安科公司(95年的)。圖3(c)(d)圖的截?cái)囝l譜數(shù)據(jù)為低頻128列。圖3(c)圖為補(bǔ)零法重建，圖3(d)圖為奇異譜分析法重建，截?cái)囝l譜數(shù)據(jù)從深圳安科公司的一幅256列頻譜數(shù)據(jù)上截下。圖象尺寸都為256×256。由實(shí)際截?cái)囝l譜數(shù)據(jù)成象情況，由圖3(a)和(b)、(c)和(d)比較得，奇異譜分析截?cái)囝l譜信號(hào)重建方法能完全消除截?cái)鄠斡?。綜上所述，奇異譜分析截?cái)囝l譜信號(hào)重建方法的重建的圖象比補(bǔ)零法重建的圖象更好是因?yàn)檠a(bǔ)零法對(duì)未采集的頻譜直接用零代替，盡管高頻部分圖象信息少、能量小，但畢竟還有小量圖象信息隱含其中，尤其高頻信息損失嚴(yán)重。毫無疑問，補(bǔ)零法完全放棄了這部分圖象信息。而本發(fā)明方法卻不一樣，它根據(jù)截?cái)囝l譜對(duì)應(yīng)圖象的邊緣特性，對(duì)那些未采集的頻譜數(shù)據(jù)進(jìn)行恢復(fù)。這就補(bǔ)救了部分本屬丟失的圖象信息和能量，很大程度上挽回了由于頻譜不足帶來的對(duì)圖象質(zhì)量的影響。尤其在奇異點(diǎn)個(gè)數(shù)較少的情況下，重建圖象的誤差僅由計(jì)算工具精度決定，即理論上講高頻分量可以準(zhǔn)確恢復(fù)。截?cái)囝l譜數(shù)據(jù)的補(bǔ)零法圖象重建是以圖象質(zhì)量(如偽影，圖象變粗糙)為代價(jià)換取成象時(shí)間縮短。本發(fā)明方法，不論含噪聲的還是無噪聲的、對(duì)實(shí)際的還是對(duì)計(jì)算機(jī)模擬的截?cái)囝l譜數(shù)據(jù)的實(shí)驗(yàn)結(jié)果，都表明了它是一種高精度的高頻信號(hào)重建方法，保證能消除重建的圖象的截?cái)鄠斡?，圖象的質(zhì)量大大優(yōu)于用現(xiàn)有方法重建的圖象。它適用于磁共振截?cái)囝l譜圖象重建和其它截?cái)囝l譜信號(hào)重建，效果顯著。表1.計(jì)算機(jī)模擬無噪聲截?cái)囝l譜數(shù)據(jù)成象的誤差、重建時(shí)間比較表</tables>表2.含5％高斯噪聲計(jì)算機(jī)模擬截?cái)囝l譜數(shù)據(jù)成象的誤差、重建時(shí)間比較表</tables>權(quán)利要求1.奇異譜分析截?cái)囝l譜信號(hào)重建方法其特征在于采用下列操作步驟為(其中，F(xiàn)=[f^(0),f^(1),...,f^(n-1)]T]]>為截?cái)囝l譜，對(duì)應(yīng)的數(shù)字信號(hào)為f＝[f(0)，f(1)，...，f(N-1)]T，n＜N，)(1)運(yùn)用下式計(jì)算各尺度下的小波系數(shù)cj,k=1NΣu=0n-1f^(u)ψ^j,k(u)‾,]]>k＝0，1，...，N-1，其中尺度j取0，1，2，3；(2)檢索極值將隨著尺度的增加而增加的極點(diǎn)，即確定奇異點(diǎn)bi，i＝1，2，...，Q。(3)用偽逆矩陣法求解以下方程組，求加權(quán)系數(shù)y＝{yb1，yb2，...，ybQ}T(4)用下式計(jì)算截?cái)囝l譜的各高頻分量f^(u)=Σi=1QybiWbi^(u),]]>u＝n、n+1，...，N-1；(5)對(duì)頻譜矢量F=[f^(0),f^(1),...,f^(n-1),f^(n),...,f^(N-1)]T]]>進(jìn)行Fourier反變換，得重建信號(hào)。2.權(quán)利要求1的奇異譜分析截?cái)囝l譜信號(hào)重建方法用于磁共振截?cái)囝l譜成象其特征是依次包括下列步驟(設(shè)磁共振截?cái)囝l譜對(duì)應(yīng)的圖象尺寸為N×N，截?cái)囝l譜的尺寸為n×N，其中n＜N)(1)從第1行到第n行頻譜數(shù)據(jù)進(jìn)行Fourier反變換(2)從第1列到第N列截?cái)囝l譜數(shù)據(jù)F=[f^(0),f^(1),...,f^(n-1)]T]]>反復(fù)執(zhí)行下例步驟1)運(yùn)用權(quán)利要求1所述操作步驟，重建截?cái)囝l譜F的信號(hào)f＝[f0，f1，...fN-1]T(圖象列矢量)2)求信號(hào)的模(圖象列矢量的各分量是個(gè)復(fù)數(shù)，取模為圖象灰度)；3)輸出圖象。全文摘要本方法關(guān)鍵是建立信號(hào)譜奇異譜分析數(shù)學(xué)模型,從采集到的已有的截?cái)囝l譜數(shù)據(jù)中抽取模型參數(shù)(奇異點(diǎn)和奇異譜分析的加權(quán)系數(shù)),再根據(jù)奇異譜分析數(shù)學(xué)模型,推出截?cái)囝l譜所隱含的高頻分量,解決截?cái)囝l譜信號(hào)重建中的偽影問題。不論含噪聲的還是無噪聲的、對(duì)實(shí)際的還是對(duì)計(jì)算機(jī)模擬的截?cái)囝l譜數(shù)據(jù)的實(shí)驗(yàn)結(jié)果,都表明了它是一種高精度的高頻信號(hào)重建方法,重建精度和速度大大優(yōu)于現(xiàn)有重建方法。它適用于磁共振截?cái)囝l譜圖象重建和其它截?cái)囝l譜信號(hào)重建。文檔編號(hào)G01N24/00GK1198530SQ98107580公開日1998年11月11日申請(qǐng)日期1998年4月20日優(yōu)先權(quán)日1998年4月20日發(fā)明者駱建華,莊天戈申請(qǐng)人:駱建華