一種超聲ct的成像方法
【專利摘要】本發明公開了一種超聲CT的成像方法。包括將成像區域網格化為M個網格點；環形探頭以成像區域為中心360°旋轉，獲得Q*N*K組數據；然后獲得該數據對應的渡越時間向量T0以及直線路徑向量L0；最后迭代計算，獲得所述M個網格點對應的速度分布V，并將所述速度分布V歸一化，完成成像。本發明通過優化環形探頭中發射陣元以及接收陣元的采集次數和采集角度，并對獲得的數據進行處理與圖像重建，從而使重建出來的速度分布更精確，以提高成像的信噪比，提高系統的成像質量。
【專利說明】
一種超聲CT的成像方法
技術領域
[0001 ]本發明屬于超聲斷層成像領域，更具體地，涉及一種超聲CT的成像方法。
【背景技術】
[0002] 超聲CT成像，即超聲斷層成像，具有超聲非侵入、無輻射和成本低的優勢，以斷層 成像的方式進行掃描，在乳腺癌等疾病早期診斷方面具有廣闊的應用前景。
[0003] 超聲CT成像分為反射式成像和透射式成像，其中，采用聲速的透射式成像是通過 重建乳腺內組織的聲速分布和衰減系數進行成像。超聲CT聲速成像利用首達波到達時間， 重建乳腺內組織的聲速分布，從而正確地區分良性腫瘤與癌。因為和正常組織與良性腫瘤 相比，癌具有較高的聲速。
[0004] 目前由在乳腺癌早期檢測領域走在前列的美國Karmano s癌癥中心研制的環形陣 列的超聲CT成像掃描器CURE可以得到亞毫米級高分辨率的圖像。這種掃描器每次使用單陣 元向環形陣列的中心發射超聲波，所有陣元接收，直到每個陣元循環發射一次，所以可以接 收到多個角度的超聲信號，得到斑點噪聲少的高分辨率圖像。雖然這種掃描器實現了對成 像物體的360度掃描，但由于環形探頭之間存在間隙，從而無法獲取到更多的成像目標的信 息，進而聲速重建時由于方程組的解并非唯一，所以求解出來的聲速分布跟實際的聲速分 布差別很大，重建圖像的質量不高。
【發明內容】
[0005] 針對現有技術的以上缺陷或改進需求，本發明提供了一種超聲CT的成像方法，其 目的在于優化環形探頭中發射陣元以及接收陣元的采集次數和采集角度，并對獲得的數據 進行處理與圖像重建，從而使重建出來的速度分布更精確，以提高成像的信噪比，提高了系 統的成像質量。
[0006] 為實現上述目的，按照本發明的一個方面，提供了一種超聲CT的成像方法，包括以 下步驟：
[0007] (1)對成像區域進行網格化處理后獲得Μ個網格點；
[0008] (2)環形探頭以成像區域為中心360°旋轉，獲得所述Μ個網格點對應的Q*N*K組數 據;其中，0 ，表示采集次數，N為環形探頭的發射陣元的數量，每個發射陣元 對應的接收陣元的數量，[.]表示向下取整；
[0009] (3)獲得所述Q*N*K組數據對應的渡越時間向量To以及直線路徑向量L〇;
[0010] (4)選取初始慢速SQ，以To為初始渡越時間、L〇為初始傳播路徑，迭代求解方程組 Lf-wdSf = dTf，dTf = Tf - Tf-i，Sf = Sf-ddSf，直至dTf <渡越時間擾動變化閾值ξ，獲得實際傳 播路徑Lf;
[0011] (5)獲得所述Μ個網格點對應的速度分布V，并將所述速度分布V歸一化，完成成像； 其中，速度分布V = 1 /s，實際慢速s = TQ/Lf。
[0012]優選地，在所述步驟⑴中，所述網格點的邊長為εχ~ε ;其中，ε ε2)，￡丄為 超聲波的半波長，ε2為成像區域中病變組織尺寸的1/4,以保證成像速度的同時，保持成像 區域中的病變組織的清晰度。
[0013] 優選地，在所述步驟（2)中，所述環形探頭旋轉的角度為#?，α為與Ν互質的整 數。
[0014] 優選地，在所述步驟(2)中，Κ為3~(Ν/2+1)。
[0015] 優選地，在所述步驟(4)中，所述初始慢速為水中的聲速。
[0016] 優選地，在所述步驟(4)中，所述渡越時間擾動變化閾值ξ為eT2~eT4。
[0017] 優選地，所述步驟(4)具體包括以下子步驟：
[0018] (4-1)令迭代次數f = 1，初始慢速Sf-i = So，初始渡越時間Tf-i = To，初始傳播路徑 Lf-i - Lo ；
[0019] (4-2)根據方程組 Lf-AdSFdThSFSf-i+dSndTFTf -Tf-!，獲得迭代渡越時間， 迭代慢速，迭代傳播路徑Lf以及傳播路徑渡越時間擾動變化dTf;
[0020] (4-3)判斷所述渡越時間擾動變化dTf是否小于渡越時間擾動變化閾值ξ，是則獲 得實際傳播路徑為迭代傳播路徑Lf，否則f = f+l，返回步驟(4-2)。
[0021] 優選地，在所述步驟(5)中歸一化的方法戈
，其中，G(i) 表示歸一化后的速度分布，V(i)表示網格點i的速度分布，i為1~M，min(V(i))為V(i)的最 小值，max(V(i))為V(i)的最大值。
[0022] 優選地，在所述步驟(5)后還包括:將所述速度分布V歸一化，并將歸一化后的速度 分布映射為灰度值。
[0023] 總體而言，通過本發明所構思的以上技術方案與現有技術相比，具有下列有益效 果：
[0024] 1、由于本發明的成像方法，使得環形探頭從多角度采集超聲數據，從面能比傳統 超聲獲取更加全面的成像數據，彌補了陣元間存在空隙而無法全面掃描到成像區域的問 題，有效緩解了聲速重建中逆問題的不適定性，從而提高了信噪比和速度重建的精度。
[0025] 2、由于環形探頭多次旋轉掃描的次數過多，會造成信息的冗余，重建時求解逆問 題的復雜度大大提高，重建效率較低等問題;而掃描次數過少，則會造成成像不清晰;本發 明根據求解慢速時的不適定性，獲得了最優的采集次數，在提高重建精度的前提下，最大限 度地保持了數據的采集效率。
【附圖說明】
[0026] 圖1是實施例1的環形陣列示意圖；
[0027] 圖2是實施例1的基于環形陣列多次旋轉掃描模式的超聲CT聲速重建流程圖；
[0028] 圖3a是實施例1的環形陣列多次旋轉掃描前重建的超聲CT聲速圖像；
[0029] 圖3b是實施例1的環形陣列旋轉一次掃描后重建的超聲CT聲速圖像。
【具體實施方式】
[0030] 為了使本發明的目的、技術方案及優點更加清楚明白，以下結合附圖及實施例，對 本發明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本發明，并 不用于限定本發明。此外，下面所描述的本發明各個實施方式中所涉及到的技術特征只要 彼此之間未構成沖突就可以相互組合。
[0031] 本發明提供了一種超聲CT的成像方法，包括以下步驟：
[0032] (1)將成像區域網格化為Μ個網格點；網格點的邊長越小，成像精度越高，然而會影 響其計算效率，因此網格點的邊長最小為超聲波的半波長，最大不超過超聲波的半波長和 成像區域中病變組織尺寸的1 /4的最大值；
[0033] (2)環形探頭以成像區域為中心360°旋轉，獲得所述Μ個網格點對應的Q*N*K組數 據;其中
，表示采集次數，N為環形探頭的發射陣元的數量，K為3~(N/2+1)，為 每個發射陣元對應的接收陣元的數量，通常接收陣元優選為以發射陣元的對面陣元為中心 呈V4~31多個陣元，以在保證成像精度的同時，盡可能的減少計算量，[.]表示向下取整;在 實踐中，通常令環形探頭以固定角度#?旋轉Q - 1次，α為與N互質的整數，為操作簡易考 慮，可令α = 1;
[0034] (3)利用加權AIC(Akaike information criterion)渡越時間提取法、采用初至波 方法(the first break method)或AIC和歸一化相結合法獲得所述Q*N*K組數據對應的渡 越時間向量Το以及直線路徑向量Lo;
[0035] (4-1)令迭代次數f = 1，初始慢速Sf-1 = So，初始渡越時間Tf-ι = Το，初始傳播路徑 LhiLo;初始慢速選取任意正數即可，為減少迭代次數，通常以水中的聲速為初始慢速；
[0036] (4-2)根據方程組 Lf-AdSf = dTf，Sf = Sf-i+dSf，dTf = Tf - Tf-!，獲得迭代渡越時間， 迭代慢速，迭代傳播路徑Lf以及傳播路徑渡越時間擾動變化dTf;
[0037] (4-3)判斷所述渡越時間擾動變化dTf是否小于渡越時間擾動變化閾值ξ，是則獲 得實際傳播路徑為迭代傳播路徑L f，否則f = f+l，返回步驟(4-2);渡越時間擾動變化閾值ξ 為eT2~eT4,要求的成像精度越高可將該閾值設置得越小；
[0038] (5)獲得所述Μ個網格點對應的速度分布V，完成成像;其中，速度分布V=l/s，實際 慢速 s = To/Lf;
[0039] (6)將所述速度分布V歸一化，并將歸一化后的速度分布映射為灰度值；歸一化可 利用公￥
其中，G(i)表不歸一化后的速度分布，V(i)表不網格點 i的速度分布，i為1~111^11(￥(1))為￥(1)的最小值，111&以￥(1))為￥(1)的最大值;映射的方 法通常采用線性映射或動態范圍壓縮法；以線性映射為例，對于一個灰度級為256的圖像， 則令max(G(i))對應映射后的灰度的最大值255，11^11(6(1))對應映射后的灰度的最小值，其 它的速度分布則按該比例映射；
[0040] 當采用線性映射時，也可以直接將歸一化前的速度分布映射為灰度值，即灰度值
[·]表示向下取整。。
[0041 ] 實施例1
[0042] 本發明提出的探頭多次旋轉掃描模式的乳腺超聲CT的聲速成像方法，通過環形探 頭作一定次數的旋轉運動，以此來獲得美國Karmanos癌癥中心的超聲CT成像掃描器CURE無 法實現的陣元間間隙所對應的成像物體的信息。同時通過規劃最優旋轉次數，降低聲速重 建時逆問題的不適定性，重建出高質量的乳腺超聲CT聲速圖像(本實施例中重建出的成像 目標的聲速與實際聲速的相對誤差減小了一倍）。
[0043] -種基于環形陣列多次旋轉掃描模式的超聲CT聲速成像方法，其步驟包括最優旋 轉掃描次數確定、環形陣列多次旋轉掃描下的數據采集、子孔徑數據選取、渡越時間向量和 初始路徑矩陣整合、折線路徑更新、聲速圖像評價標準、求逆問題，解速度分布以及圖像顯 示步驟。
[0044] (1)子孔徑數據的采集
[0045] (1-1)規劃最優旋轉掃描次數
[0046] 由于多次旋轉掃描方法使超聲環形探頭覆蓋到更多的成像區域，因此，這種掃描 方法，不僅可以采集到更豐富的信息，而且可以抑制陣元間的稀疏性造成的偽信號，從而提 高圖像的重建質量。理論上來說，只要每次旋轉的角度足夠小，旋轉掃描次數越多，采集到 的信息更全面、更豐富，但是如果掃描次數過多，會造成信息的冗余;采集時間過長，采集效 率比較低;重建時求解逆問題的復雜度大大提高，重建效率較低等問題。所以如何選取最優 的旋轉掃描次數，是本發明優先考慮的一項重要內容。本發明利用聲速重建中逆問題的不 適定性，基于最大程度地降低逆問題的不適定性來規劃最優的旋轉掃描次數。
[0047] 因此為了準確地求解成像物體的聲速分布，精確地重建圖像，我們盡可能地讓路 徑矩陣的秩等于慢度向量s的維數M，而通過多次旋轉環形探頭，可以產生更多地從發射陣 元到接收陣元的路徑，從而路徑矩陣的秩可以增大至或接近M。因此，獲得最優的旋轉掃描 次數為:……，則采集的數據的次數為Q次；
[0048] 將成像區域網格化，分成M = (m-1) * (η-1)個網格。其中，Μ為成像區域剖分的網格 數。理論上剖分的網格尺寸可以無限接近0。雖然網格尺寸越小，求解出來的聲速分布越精 確，但是網格越密，計算量越大，所以最佳的Μ由成像區域和計算效率共同決定。一般地，病 變組織的尺寸較大時，剖分不要求那么密，網格數Μ可適當降低，只要滿足網格尺寸小于病 變組織尺寸的1/4即可；當未知成像區域中病變組織尺寸時，剖分的網格大小不超過(λ/2)* (λ/2)，λ為超聲波波長。Ν為環形探頭的陣元數量，Κ為子孔徑的陣元數量;本實施例中，Ν = 72，Κ = 9，Μ = 64*64。；因為聲速重建是利用透射信號成像，子孔徑即為發射陣元對面陣元的 鄰近陣元，子孔徑接收到的信號主要包含透射信息，理論上Κ最小為3,最大為(Ν/2+1)。如果 Κ過大，不僅圖像的旁瓣水平會很高，重建的計算量也會提高;但Κ過小，能夠利用到的信息 較少，降低重建的精度。[.]表示向下取整。
[0049] (1-2)數據采集
[0050]利用Ν陣元環形探頭，循環采集成像物體與環形探頭的陣元一一對應的Ν組數據； 沿逆時針方向（或順時針方向）旋轉角度·^α，且a與Ν互質，再采集成像物體的Ν組數據;直 至沿逆時針方向采集第Q次數據，每個網格共獲得Q*N*K組子孔徑陣元接收到的信號。
[0051 ] (2)渡越時間向量和初始路徑矩陣整合：
[0052]將采集到的Q+1組子孔徑數據分別利用目前美國Karmanos癌癥中心使用的， Cuiping Li等提出的加權AIC(Akaike information criterion)渡越時間提取法(除此之 外，也可以采用初至波方法(the first break method)，Xiaolei Qu提出的AIC和歸一化互 相關結合的方法(AICNCC)等)計算每批數據的渡越時間向量，維數N*K，并按對應的旋轉掃 描的先后順序放在總的渡越時間向量矩陣To中，維數N*K*(Q+1 )，完成渡越時間向量的整 合。同時，每組子孔徑數據對應的初始路徑矩陣計算是借助成像區域網格化完成的。將成像 區域網格化，分成m行η列，這樣整個成像區域分為Μ個網格，M = (m-1) * (η-1)。由于每個網格 都對應某一個速度，而速度的倒數是慢度，所以整個成像區域的速度分布可以由一個M*1的 慢度向量s來描述。計算環陣上每個陣元在網格中的位置(第幾行第幾列）；假設超聲信號從 發射陣元沿直線路徑傳播至接收陣元，根據路徑經過的網格位置及在網格中的路徑長度表 示每個發射陣元到接收陣元的直線路徑，得到每組子孔徑數據對應的直線路徑矩陣，維數 為N*K行Μ列；最后將這Q+1個路徑矩陣也按對應的旋轉的先后順序放在總的直線路徑矩陣 Lo中，維數為N*K*(Q+1)行，Μ列，這里的直線路徑矩陣是步驟五中折線路徑更新的初始路徑 矩陣。
[0053] (3)計算速度分布V
[0054] (3-1)折線路徑更新：
[0055]超聲CT的聲速重建問題實質為求解如下一類線性方程組:Ls = T，其中，s為M*1的 慢度向量，可反映整個成像目標的速度分布，L是超聲波從發射陣元到發射陣元對面的子孔 徑的接收陣元的折線路徑矩陣，T為子孔徑中所有從發射陣元到接收陣元超聲波的渡越時 間向量。
[0056] 本發明由于規劃了最優掃描次數，從而降低了線性方程組的不適定性，從而得到 方程組更準確的解。因為聲速成像中，我們只利用到了子孔徑的信息，所以網格剖分充分密 的情況下線性方程組(1)是個欠定問題，即N*K<M。
[0057] 而由數值分析可知，設A為p階方陣，非齊次線性方程組Ax = b……（2)有唯一解等 價于系數矩陣A的秩rank(A)=p，即系數矩陣的秩等于方程組中未知數的個數。這里A對應 重建中的路徑矩陣L，X對應重建中的慢度向量s，b對應重建時的實際的渡越時間向量To，系 數矩陣A的秩對應線性無關的從發射陣元到接受陣元路徑的個數。
[0058]原理是從發射與接收點之間的假想初始路徑開始，根據最小走時準則對路徑進行 擾動，從而求出接收點處的走時及射線路徑。步驟四中計算的總的直線路徑矩陣Lo為這里 折線路徑更新的初始路徑矩陣。對由陣元i發射并且被陣元j接收的超聲波，它的傳播路徑 Lij，慢度向量Smxi和渡越時間Tij滿足:Lij(ixM)*SMxi = Tij......⑷；這樣所有的發射接收陣兀 組合的路徑矩陣L和渡越時間T滿足以下等式:L*s = T……（5);假設初始的慢度為So(任意 正數即可，但為了減少迭代次數，加快迭代速度，一般由聲波在水中的速度表示），將初始路 徑矩陣Lo和So向量相乘得到計算的渡越時間!^，然后將1^和實際的渡越時間向量作差，得出 渡越時間擾動，即cmiMo。假設在慢度分布微小擾動時，路徑矩陣沒有發生變化。解方程 組=cm，得慢度擾動dSi，又由Si=So+dSi，得更新后的慢度Si。根據慢度Si，由最短路 徑算法重新規劃每個發射陣元到接收陣元的最短路徑，然后選取每個發射陣元到子孔徑中 陣元的路徑，組成N*K行，Μ列的路徑矩陣，記為U。再將LdPSi相乘得到計算的渡越時間T 2， 計算新的渡越時間擾動，即dTsiTs-Ti，同樣利用在慢度分布微小擾動時，路徑矩陣沒有發 生變化的假設，解方程組Lo*dS 2 = dT2，得慢度擾動dS2，又由S2 = Si+dS〗，得更新后的慢度S2， 重復之前的計算，不斷更新折線路徑，不斷更新慢度分布。
[0059] (3-2):聲速圖像評價標準:該評價標準指的是步驟五中每次更新折線路徑時得到 的渡越時間擾動變化dT f的范數上限ξ，計算的渡越時間和實際的渡越時間之差稱為渡越時 間的擾動。渡越時間擾動變化的上限是一個很小的數，ξ-般小于eT 2(在本實施例中| = e _2)，當渡越時間擾動變化不超過上限時，這時就認為達到我們所要求的精度，就可以進行下 一步圖像顯示;如果渡越時間擾動變化超過規定的上限，則繼續更新折線路徑，更新f-Ι次， 直到擾動變化dT f不超過上限，這時的折線路徑矩陣記為Lf，對應的慢度為&。
[0060] (3-3)求逆問題，解速度分布:通過步驟五和步驟六中得到的最終的折線路徑矩陣 Lf和步驟四中得到的渡越時間向量Το，通過一種成熟的方法-TVAL3方法(其他還有非線性共 輒梯度方法，聯合代數重建算法等)求解逆問題:L f*S = T〇……（6)得出慢度分布s，利用V = 1/8，將8轉化成速度向量￥，維數為(111-1)*(11-1)。
[0061] (4)圖像顯示:該步驟包括將步驟七中得到的速度向量歸一化及灰度映射，最終得 到超聲CT圖像。為了提高圖像的對比度，我們采用將速度向量V歸一化。將V的最小值min (V)，最大值max(V)取出來，然后按照(7)中函數將其映射到[0，1]，并最終生成(m-l)*(n-l) 的超聲CT圖像。灰度映射采用簡單的線性映射(其他還有動態范圍壓縮法等），即成比例地 將最小的速度映射到〇,將最大的速度映射到255或511 (分別對應圖像顯示中常用的256灰 度級或512灰度級的圖像）。
[0062]
[0063]從圖3中我們可以看出旋轉掃描一次對聲速重建圖像的影響，可以看出，經一次掃 描，聲速重建后的圖像質量即獲得了極大的提升。
[0064]本領域的技術人員容易理解，以上所述僅為本發明的較佳實施例而已，并不用以 限制本發明，凡在本發明的精神和原則之內所作的任何修改、等同替換和改進等，均應包含 在本發明的保護范圍之內。
【主權項】
1. 一種超聲CT的成像方法，其特征在于，包括以下步驟： (1) 對成像區域進行網格化處理后獲得M個網格點； (2) 環形探頭以成像區域為中心360°旋轉，獲得所述M個網格點對應的Q*N*K組數據;其 中，表示采集次數，N為環形探頭的發射陣元的數量，K為每個發射陣元對應的 接收陣元的數量，[.]表示向下取整； (3) 獲得所述Q*N*K組數據對應的渡越時間向量To以及直線路徑向量L0; (4) 選取初始慢速So，以To為初始渡越時間、Lo為初始傳播路徑，迭代求解方程組Lf-AdSf = dTf，dTf = Tf - Tf-! ,Sf = Sf-WdSf，直至dTf<渡越時間擾動變化閾值ξ，獲得實際傳播路徑 Lf; (5) 獲得所述M個網格點對應的速度分布V，完成成像;其中，速度分布V= Ι/s，實際慢速 s = To/Lf 〇2. 如權利要求1所述的成像方法，其特征在于，在所述步驟（1)中，所述網格點的邊長為 ει~ε ;其中，E=Iiiax(EU2) J1為超聲波的半波長，ε2為成像區域中病變組織尺寸的1/4。3. 如權利要求1所述的成像方法，其特征在于，在所述步驟(2)中，所述環形探頭旋轉的 角度戈，α為與N互質的整數。4. 如權利要求1所述的成像方法，其特征在于，在所述步驟(2)中，K為3~(Ν/2+1)。5. 如權利要求1所述的成像方法，其特征在于，在所述步驟(4)中，所述初始慢速為水中 的聲速。6. 如權利要求1所述的成像方法，其特征在于，在所述步驟(4)中，所述渡越時間擾動變 化閾值ξ為e_2~e_4。7. 如權利要求1所述的成像方法，其特征在于，所述步驟(4)具體包括以下子步驟： (4-1)令迭代次數f = 1，初始慢速Sf-i = So，初始渡越時間Tf-i = To，初始傳播路徑Lf-i = Lo； (4-2)根據方程組Lf-AiSf = dTf ,Sf = Sf-WdSf，dTf = Tf- Τη，獲得迭代渡越時間，迭代 慢速，迭代傳播路徑Lf以及傳播路徑渡越時間擾動變化dTf; (4-3)判斷所述渡越時間擾動變化dTf是否小于渡越時間擾動變化閾值ξ，是則獲得實際 傳播路徑為迭代傳播路徑Lf，否則f = f+l，返回步驟(4-2)。8. 如權利要求1所述的成像方法，其特征在于，在所述步驟（5)中歸一化的方法為其中，G(i)表不歸一化后的速度分布，V(i)表不網格點i的速度 分布，i為1~M，min(V(i))為V(i)的最小值，max(V(i))為V(i)的最大值。9. 如權利要求1-8中任意一項所述的成像方法，其特征在于，在所述步驟(5)后還包括： 將所述速度分布V歸一化，并將歸一化后的速度分布映射為灰度值。
【文檔編號】A61B8/15GK105997153SQ201610559079
【公開日】2016年10月12日
【申請日】2016年7月15日
【發明人】尉遲明, 丁明躍, 婁翠娟, 王珊珊, 宋俊杰, 李春雨, 方小悅, 彭楊
【申請人】華中科技大學