本發(fā)明涉及一種MRI成像技術應用中非線性最小二乘算法的優(yōu)化改進方法，特別涉及一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法。
背景技術：
：基于多加權b值的擴散磁共振成像技術越來越成為擴散磁共振技術的發(fā)展趨勢，包括擴散峰度成像、擴散譜成像等。同時，模型的階數也在不斷提高，例如從二階的擴散張量成像到四階的擴散峰度成像模型。因而模型的對于數據信噪比的敏感性也越來越敏感，因而在實際應用中，由于模型擬合的誤差出現(xiàn)奇異值、擬合不穩(wěn)定、無意義解的情況也逐漸增多，影響技術的穩(wěn)定發(fā)揮和應用。本技術發(fā)明旨在通過一定的算法改進優(yōu)化，提高多b值模型的尋優(yōu)求解穩(wěn)定性，進一步促進MRI成像技術的發(fā)展和應用。技術實現(xiàn)要素：針對現(xiàn)有技術存在的技術問題，本發(fā)明提供一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法，該方法是針對多b值擴散加權構建擬合模型，通過最小二乘算法的優(yōu)化提高模型尋優(yōu)求解的穩(wěn)定性，進而使得MRI成像中參數圖像的估計更加穩(wěn)定，降低奇異值或無意義值的出現(xiàn)，進一步促進MRI成像技術的發(fā)展和應用。為了解決現(xiàn)有技術中存在技術問題，本發(fā)明采用如下步驟予以實施：一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法，包括下列步驟：步驟一，通過多b值加權數據建立擬合模型；步驟二，通過多b值加權數據估算信噪比；步驟三，根據信噪比計算殘差修正權值αi，步驟四，建立擬合模型初始值；步驟五，根據殘差修正權值αi計算擬合模型殘差修正值；步驟六，如果擬合模型殘差修正值符合收斂條件，則通過擬合模型的最優(yōu)解計算提高其穩(wěn)定性；否則，返回步驟四繼續(xù)尋求擬合模型的最優(yōu)解。所述步驟一中的擬合模型建立Si＝f(bi，β)，，其中β為待求解的未知參數，Si>0，i＝0,1,2,3…；b0＝0；bi>0，i＝1,2,3…。所述步驟二中信噪比的估算是包括全局信噪比和局部信噪比，即SNRi和SNR0。所述步驟三中殘差修正權值αi采用如下公式：這里，b0＝0；bi>0,i＝1,2,3,….,bi常取值500,1000,1500,2000,2500,…。所述步驟五中擬合模型殘差修正值計算采用如下公式：E(β^)=Σi=0nαi×(Si-f(bi,β^))2.]]>所述步驟六中擬合模型最優(yōu)解計算采用如下公式：min.E(β^)=Σi=0nαi×(Si-f(bi,β^))2.]]>本發(fā)明有益效果：第一，本發(fā)明通過估計不同加權b值數據的信噪比加權，削減或者消除不同加權b值信噪比差異，進而減弱了噪聲和擬合誤差對模型擬合穩(wěn)定性的影響。第二，如圖2所示，本發(fā)明通過在模型優(yōu)化求解過程中引入對殘差平方和進行修正的權值向量，使得曲線擬合的參數結果更依賴于高信噪比的加權b值，即有限保證更高信噪比(較低)的b值點數據，從而減小高b值中噪聲對曲線擬合的影響，使得參數的求解更加穩(wěn)定。第三，如圖3所示，奇異值的情況多出現(xiàn)在灰質和灰質區(qū)域，本發(fā)明將本來較高的灰質白質峰度值估計成較低值甚至零值。這里分別給出本發(fā)明方法對灰質(圖4)、白質(圖5)奇異值以及正常(圖6)組織的效果。圖4、圖5和圖6中的E和Ec，分別是修正前的求解函數曲面(式5)和修正后的求解函數曲面(式6)，圖中的白色“+”是兩種方法得到的解的位置；相應地，圖中下方給出了不同方法得到的擬合曲線，其中藍色均為使用本發(fā)明方法得到的曲線擬合結果。此三圖結果明顯展示了本發(fā)明方法的效果：能夠很好解決擬合過程中出現(xiàn)的奇異值現(xiàn)象，同時不影響原本正確穩(wěn)定的擬合結果。第四，如圖7所示，給出了現(xiàn)有方法和本發(fā)明方法得到的擴散峰度成像參數圖的結果，其中上方兩圖為現(xiàn)有方法得到的參數圖，下方的兩圖為使用本發(fā)明方法得到的參數圖。附圖說明圖1是本發(fā)明一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法流程圖。圖2是本發(fā)明一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法對殘差以及擬合結果的影響圖。圖3是本發(fā)明一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法的模型擬合結果極易出現(xiàn)的奇異值示意(黑點)圖。圖4是本發(fā)明一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法的灰質或腦脊液的奇異值圖。圖5本發(fā)明一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法的白質部位奇異值圖。圖6是本發(fā)明一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法的優(yōu)化方案對正常部位的圖。圖7是本發(fā)明一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法的優(yōu)化方案的實際效果展示圖(左：平均峰度；右：徑向峰度)。具體實施方式下面結合附圖對本發(fā)明做出詳細地說明：如圖1所示，本發(fā)明提供一種基于信噪比加權的多b值擴散磁共振成像優(yōu)化方法，包括下列步驟：步驟一101，通過多b值加權數據建立擬合模型；所述步驟一中的擬合模型建立Si＝f(bi，β)，，其中β為待求解的未知參數，Si>0，bi>0，i＝0,1,2,3…。在實際中，多b值加權成像在目前的MRI技術應用中已經較為常見，其數據包括一個b＝0的無擴散加權的基準信號(圖像)，和若干非零b值下的擴散加權成像信號(圖像)，且這些非零b值擴散加權信號(圖像)來自同一個對象。這里設基準信號為強度為S0＝S(b0)，而加權b值的信號為Si＝S(bi)(i＝1,2,3，…)。對于某個對象(仿體實物、動植物等所有非金屬物)進行多b值擴散加權成像時，信號S的強度隨著b值變大而不斷減小，構成的曲線稱之為衰減曲線，實際應用中用不同的模型進行擬合并定量描述對象的結構特性。這里將所有模型一般抽象地用Si＝f(bi，β)表示，其中β為待求解的未知參數，Si>0，bi>0，i＝0,1,2,3…。步驟二102，通過多b值加權數據估算信噪比；所述步驟二中信噪比的估算是包括全局信噪比和局部信噪比，即SNRi和SNR0。這里需要說明，在相同的機器環(huán)境中，不同b值數據(信號或者圖像)中噪聲的水平是相同，而由于不同b值加權得到的信號強度不同，使得不同b值加權數據的信噪比不同。這里需要首先估計不同b值加權圖像的信噪比，包括全局信噪比和局部信噪比。其中，全局信噪比SNRi＝MeanSignal/MeanNoise。這里MeanSignal為第i個b值圖像中信號(檢測對象或目標區(qū)域的信號)的均值，MeanNoise為第i個b值圖像中空氣噪聲信號(非檢測對象或非目標區(qū)域的信號)的均值。局部信噪比SNRi通過局部空間(或者圖像的局部區(qū)域/鄰域)信號的均值除以信號的標準差來估計。步驟三103，根據信噪比計算殘差修正權值αi，所述步驟三中殘差修正權值αi采用如下公式：這里，b0＝0；bi＞0.i＝1，2，3，…，bi常取值500，1000，1500，2000，2500，··…（1）在數據曲線的擬合精度一定程度上依賴于采樣點(已知點)數據的可靠性。因為不同加權b值下得到的信號其信噪比不同，因而其中有用信號或者數據的可信程度也不一樣，這會對曲線或這模型的擬合造成一定的偏差。具體來說，采樣點不等同的信噪比使得噪聲對曲線擬合精度的影響更顯著了，要想控制這一影響，首先需要根據信噪比估計不同b值下信號的可靠程度。步驟四104，建立擬合模型初始值；是對擬合模型中β的某個具體值(初始值)，采用如下公式計算：Si=f(bi,β^)+μ---(2)]]>其殘差平方和為：E(β^)=Σi=0n(Si-f(bi,β^))2---(3)]]>步驟五105，根據殘差修正權值αi計算擬合模型殘差修正值；所述步驟五中擬合模型殘差修正值計算采用如下公式：E(β^)=Σi=0nαi×(Si-f(bi,β^))2.---(4)]]>步驟六(106，107)，如果擬合模型殘差修正值符合收斂條件，則通過擬合模型的最優(yōu)解計算提高其穩(wěn)定性；否則，返回步驟四繼續(xù)尋求擬合模型的最優(yōu)解。所述步驟六中擬合模型最優(yōu)解計算采用如下公式：minE(β^)=Σi=0nαi×(Si-f(bi,β^))2.---(5)]]>上式子取極小值的一階條件：dE(β^)dβ^=-2Σi=1nαi×(Si-f(bi,β^))×(-df(bi,β^)dβ^)=0---(6)]]>引入權值前后對于尋找最優(yōu)解的過程和方法沒有任何影響，只是殘差平方和發(fā)生了改變，形象來說是最優(yōu)解的位置發(fā)生了變化。改進的殘差平方和通過對每個點等權值的累加，而改進后的殘差平方和為不同點估計誤差的加權平方和，因此在此新的殘差和函數上尋求最優(yōu)解。特別地，該殘差和函數的維度與β的位置參數的個數一致。若β有兩個未知參數，殘差和函數的緯度為2，為一個二次曲面；而殘差函數的極小值一階條件也將變成兩個偏一階導聯(lián)立等于零。本發(fā)明中擬合模型尋求最優(yōu)解返回步驟四的過程是通過迭代尋優(yōu)方法，其中以β＝(β0，β1)為例，取極小值的一階條件為：dE(β0,β1)dβ0=-2Σi=1nαi×(Si-f(bi,β0,β1))×(-df(bi,β0,β1)dβ1)=0---(7)]]>dE(β0,β1)dβ1=-2Σi=1nαi×(Si-f(bi,β0,β1))×(-df(bi,β0,β1)dβ1)=0---(8)]]>擴散磁共振技術中多b值的擴散衰減模型多為非線性的，如二次曲線或者雙指數曲線等，因而基本采用非線性最小二乘的求解思路。在具體實踐中，也會加入一些對未知參數的約束條件，而采用一些不同的求解過程，大都也是非線性最小二乘的變式。非線性最小二乘法的思路是，通過泰勒級數將均值函數展開為線性模型。即，只包括一階展開式，高階展開式都歸入誤差項。然后再進行OLS回歸，將得到的估計量作為新的展開點，在對線性部分進行估計。如此往復，直至收斂。這里給出高斯-牛頓(Gauss-Newton)迭代法求解非線性最小二乘問題的數值求解過程。對原始模型展開泰勒級數，取一階近似值：f(bi,β^)≈f(bi,β^(0))+df(bi,β^)dβ^|β^(0)(β^-β^(0))---(9)]]>zi(β^)=df(bi,β^)dβ^---(10)]]>E(β^)=Σi=0nαi×(Si-f(bi,β^)-zi(β^(0))(β^-β^(0)))2=Σi=1nαi×(Si-f(bi,β^)+zi(β^(0))β^(0)-zi(β^(0))β^)2=Σi=1nαi×(St~(β^(0))-zi(β^(0))β^)2---(11)]]>構造并估計線性偽模型：St~(β^(0))=zi(β^(0))β^+ϵi---(12)]]>Et~(β^(1))=Σi=0nαi×(St~(β^(0))-zi(β^(0))β^(1))2---(13)]]>估計得到第一次迭代值并繼續(xù)迭代。完成迭代過程如下：第一步：給出參數估計值的初始值將在處展開泰勒級數，取一階近似值；第二步：計算和的樣本觀測值；第三步：采用普通最小二乘法估計模型得到β的估計值第四步：用代替第一步中的重復這些過程，直至收斂。收斂條件的設定有較大的自由度，例如設置為100次連續(xù)迭代殘差等于0或者不變可認為收斂。本發(fā)明方法是力圖解決擴散磁共振成像多b值模型擬合問題中，不同b值信噪比不同尤其是高b值信噪比很低對模型擬合結果精度和穩(wěn)定性造成破壞的問題。以實際應用中擴散峰度成像(diffusionkurtosisimaging，DKI)模型為例，其擬合函數為：f(bi，β)＝-bi×D+(1/6)×(bi×D)2×K＝-bi×D+(1/6)×bi2×D×K這里β＝[DK]。因此該求解問題是一個2參數的函數(式5、6)，其問題為在一個二維曲面上求取最小值，如圖3、4、5中的E和Ec。由于K是一個四階統(tǒng)計量，對模型中的噪聲和誤差異常敏感，因而極易出現(xiàn)如圖2所示的奇異值點。圖4、5中，可以看出求解函數的解所在范圍在K的維度上跨度更大，更體現(xiàn)了擬合誤差在K值估計中的不穩(wěn)定性。而本發(fā)明方法在DKI應用中的目的就是解決DKI對峰度估計的不穩(wěn)定性問題。DKI使用應用中至少需要不少于2個的非零b值(b值量級為103，單位為s/mm2)，而且其中必須有不低于b＝2000的高b值，這里，高b值必將引入較低信噪比的擬合點，影響峰度的估計。采用本發(fā)明方法能夠明顯改善提高DKI模型的峰度估計，如圖6所示。上述實例僅用于說明本發(fā)明，其中各部件的結構、材料、連接方式都是可以有所變化的，凡是在本發(fā)明技術基礎上進行的等同變換和改進，均不應該排除在本發(fā)明的保護范圍之外。當前第1頁1 2 3