本發(fā)明屬于核磁共振領域，更具體地，涉及一種基于相位編碼的梯度勻場方法。

背景技術：
核磁共振(NuclearMagneticResonance)作為一種重要的分析檢測工具，已經在生物大分子、代謝組學、醫(yī)藥研發(fā)等領域，取得了廣泛的應用。隨著核磁共振研究領域的深入發(fā)展，為了提高分子更為精細的結構，對核磁譜圖的分辨率也要求越來越高，而磁場的均勻性是影響譜圖分辨率的重要因素。因此，如何提高磁場的均勻性是獲得高分辨核磁共振譜圖的重要課題。梯度勻場是一種借助脈沖梯度場采集場圖矩陣和當前磁場分布數據，然后通過勻場線圈對磁場不均勻性進行補償的方法，是核磁共振應用較為廣泛的勻場方法。核磁科研工作者們已經在基于梯度勻場的方法上進行了一系列卓有成效的工作，vanZijl等提出了核磁共振的三維梯度勻場方法(Vanzijl,P.C.M.,etal.,OptimizedShimmingforHigh-ResolutionNMRUsing3-DimensionalImage-BasedField-Mapping.JournalofMagneticResonanceSeriesA,1994.111(2):p.203-207)，Weiger提出了一種基于譜圖線性優(yōu)化的梯度勻場方法(Weiger,M.,T.Speck,andM.Fey,GradientShimmingwithSpectrumOptimisation.JournalofMagneticResonance,2006.182(1):p.38-48)，Nishiyama等提出一種基于魔角旋轉的梯度勻場方法(Nishiyama,Y.,Y.Tsutsumi,andH.Utsumi,MagicShimming:GradientShimmingWithMagicAngleSampleSpinning.JournalofMagneticResonance,2012.216:p.197-200)。梯度勻場雖然已經在核磁共振的勻場中取得了重要的應用，但是在基礎場均勻性很差的情況下勻場效果較差。原因主要有兩個：1.基礎場均勻性很差時，梯度場的線性度會產生一定的偏差，使得梯度回波得到的圖像形成一定程度的扭曲；2.傳統梯度勻場存在準備梯度和讀梯度，需要先散相后重聚，編碼需要的時間較長，而較差的磁場均勻性使得脈沖激發(fā)出來的核磁信號衰減速度很快，使得時域回波信號的信噪比大為降低。

技術實現要素：
針對現有技術的以上缺陷或改進需求，本發(fā)明提供了一種基于相位編碼的梯度勻場方法，由于其利用相位編碼不受磁場不均勻性的影響，且編碼所需時間短的優(yōu)勢，減小由較差的磁場均勻性所帶來的圖像扭曲和信號衰減，可得到更為準確、信噪比更高的場圖矩陣和當前磁場分布數據，從而達到更好的勻場效果，由此解決現有技術對于基礎場均勻性較差的情況下勻場效果不好的技術問題。為實現上述目的，按照本發(fā)明的一個方面，提供了一種基于相位編碼的梯度勻場方法，其特征在于，包括以下步驟：(1)獲取場圖矩陣；(2)采集并計算當前磁場分布數據；(3)根據所述場圖矩陣和當前場分布數據，采用Tikhonov規(guī)則化的奇異值分解法，得到各勻場線圈的電流改變量，并根據所述各勻場線圈的電流改變量重新配置勻場線圈電流；(4)根據步驟(3)中獲得的各勻場線圈的電流改變量，計算相位殘差：當相位殘差滿足收斂條件時結束，否則將各勻場線圈的電流修改后的磁場作為當前場，重復步驟(2)至(4)。優(yōu)選地，所述梯度勻場方法，其步驟(1)包括以下子步驟：(1-1)采樣：核磁共振儀器在基礎場條件下，載入場圖單點成像脈沖序列，并按照預設的場圖第一及第二采樣參數進行采樣，得到第一基礎場時域回波信號EchoAbase及第二基礎場時域回波信號EchoBbase；(1-2)修改：依次修改通過各組勻場線圈的勻場電流，每修改一組電流值，按照(1-1)中載入的脈沖序列和采樣參數進行采樣，得到第一場圖時域回波信號組EchoA0shimchannel及第二場圖時域回波信號組EchoB0shimchannel；(1-3)計算場圖矩陣FieldMap(shimchannel)：其中，Δφbase為基礎磁場的相位差數據，Δφshimchannel為各勻場線圈電流相對于基礎場修改后對應的相位差數據。優(yōu)選地，所述梯度勻場方法，其步驟(1-1)其所述場圖第一及第二采樣參數包括：場圖第一及第二相位編碼時間tp01及tp02、場圖第一及第二梯度強度數組gp01及gp02。所述場圖第一及第二采樣參數符合以下條件：其中，i表示場圖相位編碼的步數，ΔG01和ΔG02表示場圖相位編碼對應梯度的步進，且滿足ΔG01*tp01＝ΔG02*tp02。優(yōu)選地，所述梯度勻場方法，其基礎磁場的相位差數據Δφbase，按照如下方法計算：Δφbase＝φ(r,tp02)base-φ(r,tp01)base其中，φ(r,tp01)base為基礎場第一相位數據，由EchoAbase經快速傅里葉變換后計算相位得到；φ(r,tp02)base為基礎場第二相位數據，由EchoBbase經快速傅里葉變換后計算相位得到。優(yōu)選地，所述梯度勻場方法，其各勻場線圈電流相對于基礎場修改后對應的相位差數據Δφshimchannel，按照如下方法計算：Δφshimchannel＝φ(r,tp02)shimchannel-φ(r,tp01)shimchannel，其中，φ(r,tp01)shimchannel為場圖第一采樣參數下勻場線圈電流修正后的相位數據，由EchoA0shimchannel經快速傅里葉變換后計算相位得到；為場圖第二采樣參數下勻場線圈電流修正后的相位數據，由EchoB0shimchannel經快速傅里葉變換后計算相位得到。優(yōu)選地，所述梯度勻場方法，其步驟(2)包括以下子步驟：(2-1)采樣：核磁共振儀器載入單點成像序列按照預設的第一及第二采樣參數進行采樣，得到第一及第二當前場時域回波信號EchoAcurr及EchoBcurr；(2-2)計算當前場分布數據：所述當前場分布數據B0curr按照如下方法計算：其中為當前場的相位差數據，按照如下方法計算：其中，φ(r,tp1)curr為當前場第一相位數據，由EchoAcurr經快速傅里葉變換后計算相位得到；φ(r,tp2)curr為當前場第二相位數據，由EchoBcurr經快速傅里葉變換后計算相位得到。優(yōu)選地，所述梯度勻場方法，其第一及第二采樣參數包括：第一及第二相位編碼時間tp1及tp2、第一及第二梯度強度數組gp1及gp2；所述第一及第二采樣參數符合以下條件：其中，n表示當前場相位編碼的步數，ΔG1和ΔG2表示當前場相位編碼對應梯度的步進，且滿足ΔG1*tp1＝ΔG2*tp2。優(yōu)選地，所述梯度勻場方法，其步驟(3)具體為：采用Tikhonov規(guī)則化的奇異值分解法解如下方程組：FieldMap*x＝B0curr其中x為勻場線圈電流的改變量。優(yōu)選地，所述梯度勻場方法，其收斂條件為：相位殘差的L2范數小于預設的L2范數閾值。總體而言，通過本發(fā)明所構思的以上技術方案與現有技術相比，能夠取得下列有益效果：實驗證實，本發(fā)明與傳統的梯度回波勻場方法相比，由于采用基于單點成像(Single-PointImaging)的純相位編碼采樣方法，能明顯的提高勻場效果。附圖說明圖1為本發(fā)明方法的總流程圖。圖2為本發(fā)明實施例所用脈沖序列的示意圖。圖3-1為梯度回波方法獲取的場圖矩陣。圖3-2為本發(fā)明實施例獲取的場圖矩陣。圖4-1為度回波方法獲取的梯度Profile。圖4-2為本發(fā)明實施例獲取的梯度Profile。圖5-1為初始磁場條件下水樣的1H譜線形。圖5-2為采用傳統的梯度回波自動勻場后水樣的1H譜線形。圖5-3為本發(fā)明實施例勻場后水樣的1H譜線形。具體實施方式為了使本發(fā)明的目的、技術方案及優(yōu)點更加清楚明白，以下結合附圖及實施例，對本發(fā)明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。此外，下面所描述的本發(fā)明各個實施方式中所涉及到的技術特征只要彼此之間未構成沖突就可以相互組合。本發(fā)明提供的基于相位編碼的梯度勻場方法，總體流程圖如圖1所示，包括以下步驟：(1)獲取場圖矩陣；優(yōu)選按照如下方法獲取場圖矩陣：(1-1)采樣：核磁共振儀器在基礎場條件下，載入場圖單點成像脈沖序列，并按照預設的場圖第一及第二采樣參數進行采樣，得到第一基礎場時域回波信號EchoAbase及第二基礎場時域回波信號EchoBbase；所述場圖第一及第二采樣參數包括：場圖第一及第二相位編碼時間tp01及tp02、場圖第一及第二梯度強度數組gp01及gp02。所述場圖第一及第二采樣參數符合以下條件：其中，i表示場圖相位編碼的步數，ΔG01和ΔG02表示場圖相位編碼對應梯度的步進，且滿足ΔG01*tp01＝ΔG02*tp02。對于按照場圖第一采樣參數采集的數據點，對應于tp01、gp01按步進變化采集的所有數據點組成的第一基礎場時域回波信號EchoAbase；對于按照場圖第二采樣參數采集的數據點，對應于tp02、gp02按步進變化采集的所有數據點組成的第二基礎場時域回波信號EchoBbase。(1-2)修改：依次修改通過各組勻場線圈的勻場電流，每修改一組電流值，按照(1-1)中載入的脈沖序列和采樣參數進行采樣，得到第一場圖時域回波信號組EchoA0shimchannel及第二場圖時域回波信號組EchoB0shimchannel。(1-3)計算場圖矩陣：場圖矩陣FieldMap(shimchannel)按照如下方法計算：其中，Δφbase為基礎磁場的相位差數據，按照如下方法計算：Δφbase＝φ(r,tp02)base-φ(r,tp01)base其中，φ(r,tp01)base為基礎場第一相位數據，由EchoAbase經快速傅里葉變換后計算相位得到；φ(r,tp02)base為基礎場第二相位數據，由EchoBbase經快速傅里葉變換后計算相位得到；Δφshimchannel為各勻場線圈電流相對于基礎場修改后對應的相位差數據，按照如下方法計算：Δφshimchannel＝φ(r,tp02)shimchannel-φ(r,tp01)shimchannel其中，φ(r,tp01)shimchannel為勻場線圈電流修正后第一采樣參數下的相位數據，由EchoA0shimchannel經快速傅里葉變換后計算相位得到；φ(r,tp2)shimchannel為勻場線圈電流修正后第二采樣參數下的相位數據，由EchoB0shimchannel經快速傅里葉變換后計算相位得到。(2)采集并計算當前磁場分布數據；(2-1)采樣：核磁共振儀器載入單點成像序列按照預設的第一及第二采樣參數進行采樣，得到第一及第二當前場時域回波信號；所述第一及第二采樣參數包括：第一及第二相位編碼時間tp1及tp2、第一及第二梯度強度數組gp1及gp2。所述當前場下的第一及第二采樣參數符合以下條件：其中，n表示當前場相位編碼的步數，ΔG1和ΔG2表示當前場相位編碼對應梯度的步進，且滿足ΔG1*tp1＝ΔG2*tp2對于按照第一采樣參數采集得到的數據點，對應于tp1、gp1按步進變化采集的所有數據點組成第一當前場時域回波信號EchoAcurr；對于按照第二采樣參數采集得到的數據點，對應于tp2、gp2按步進變化采集的所有數據點組成第二當前場時域回波信號EchoBcurr。優(yōu)選第一及第二采樣參數與場圖第一及第二采樣參數相同。(2-2)計算當前場分布數據：所述當前場分布數據B0curr按照如下方法計算：其中Δφcurr為當前場的相位差數據，按照如下方法計算：Δφcurr＝φ(r,tp2)curr-φ(r,tp1)curr其中，φ(r,tp1)curr為當前場第一相位數據，由EchoAcurr經快速傅里葉變換后計算相位得到；φ(r,tp2)curr為當前場第二相位數據，由EchoBcurr經快速傅里葉變換后計算相位得到。(3)根據所述場圖矩陣和當前場分布數據，采用Tikhonov規(guī)則化的奇異值分解法，得到各勻場線圈的電流改變量，并根據所述各勻場線圈的電流該變量重新配置勻場線圈電流。具體步驟為：采用Tikhonov規(guī)則化的奇異值(singularvaluedecomposition)分解法解如下方程組：FieldMap*x＝B0curr其中x為勻場線圈電流的改變量，將勻場線圈電流值修改x后，重新配置到勻場線圈中。(4)根據步驟(3)中獲得的各勻場線圈的電流改變量，計算磁場的相位殘差：當相位殘差滿足收斂條件時結束，否則將各勻場線圈的電流修改后的磁場作為當前場，重復步驟(2)至(4)。所述收斂條件優(yōu)選為：相位殘差的L2范數小于預設的L2范數閾值，即||B0curr-FieldMap*x||2<ε其中，|B0curr-FieldMap*x|為相位殘差，ε為預設的L2范數閾值。以下為實施例：本領域的技術人員容易理解，以上所述僅為本發(fā)明的較佳實施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內所作的任何修改、等同替換和改進等，均應包含在本發(fā)明的保護范圍之內。本實施例使用中國科學院武漢物理與數學研究所自主研制的500MHz核磁共振波譜儀，配備28組勻場線圈。下面通過實施例，結合附圖，對本發(fā)明的技術方案作進一步的詳細描述。一種基于相位編碼的梯度勻場方法，總體流程圖如圖1所示，該方法總體包括以下步驟：(1)獲取場圖矩陣；場圖矩陣的獲取方法如下：(1-1)采樣：核磁共振儀器在基礎場條件下，載入單點成像序列(Single-PointImaging序列，如圖2所示)并按照預設的場圖第一及第二采樣參數進行采樣，得到第一基礎場時域回波信號EchoAbase及第二基礎場時域回波信號EchoBbase；所述場圖第一及第二采樣參數包括：場圖第一及第二相位編碼時間tp01及tp02、場圖第一及第二梯度強度數組gp01及gp02。所述場圖第一及第二采樣參數符合以下條件：其中，i表示場圖相位編碼的步數，ΔG01和ΔG02表示場圖相位編碼對應梯度的步進，且滿足ΔG01*tp01＝ΔG02*tp02對于按照場圖第一采樣參數采集的數據點，對應于tp01、gp01按步進變化采集的所有數據點組成的第一基礎場時域回波信號EchoAbase；對于按照場圖第二采樣參數采集的數據點，對應于tp02、gp02按步進變化采集的所有數據點組成的第二基礎場時域回波信號EchoBbase。在本實施例中：采樣參數主要包括譜寬sw＝50kHz，每次掃描采樣點數fidpoints＝256，n＝64，tp1＝1ms，ΔG1＝0.130Gauss/cm，tp2＝2ms，ΔG2＝0.065Gauss/cm。(1-2)修改：依次修改通過各組勻場線圈的勻場電流，每修改一組電流值，按照(1-1)中載入的脈沖序列和采樣參數進行采樣，得到場圖第一時域回波信號組EchoA0shimchannel及場圖第二時域回波信號組EchoB0shimchannel。在本實施例中：共修改Z1～Z7共七組勻場線圈的勻場電流，各組勻場線圈電流值對應的變化量如下：(1-3)計算場圖矩陣：依據步驟(1-1)采得的數據，計算基礎場相位差數據Δφbase，計算方法如下：Δφbase＝φ(r,tp02)base-φ(r,tp01)base其中，φ(r,tp01)base為基礎場第一相位數據，由EchoAbase經快速傅里葉變換后計算相位得到；φ(r,tp02)base為基礎場第二相位數據，由EchoBbase經快速傅里葉變換后計算相位得到。場圖矩陣FieldMap(shimchannel)按照如下方法計算：其中，Δφshimchannel為各勻場線圈電流修改后對應的相位差數據，按照如下方法計算：Δφshimchannel＝φ(r,tp02)shimchannel-φ(r,tp01)shimchannel其中，φ(r,tp01)shimchannel為勻場線圈電流修正后第一采樣參數下的相位數據，由EchoA0shimchannel經快速傅里葉變換后計算相位得到；為勻場線圈電流修正后第二采樣參數下的相位數據，由EchoB0shimchannel經快速傅里葉變換后計算相位得到。在本實施例中：本方法采集得到的場圖矩陣FieldMap(shimchannel)如圖3-2所示，信噪比較圖3-1梯度回波采得的場圖矩陣有了明顯的提高。(2)測量并計算當前磁場分布數據；當前場的測量方法如下：(2-1)采樣：核磁共振儀器載入單點成像序列(Single-PointImaging序列，如圖2所示)，按照預設的第一及第二采樣參數進行采樣，得到第一當前場時域回波信號EchoAcurr及第二當前場時域回波信號EchoBcurr；所述第一及第二采樣參數包括：第一及第二相位編碼時間tp1及tp2、第一及第二梯度強度數組gp1及gp2，取值與步驟(1-1)類似。對于當前場下載入單點成像序列第一采樣參數采集的數據點，對應于tp1、gp1按步進變化采集的所有數據點組成EchoAcurr；對于當前場下載入單點成像序列第二采樣參數采集的數據點，對應于tp2、gp2按步進變化采集的所有數據點組成EchoBcurr。(2-2)計算當前場分布數據：所述當前場分布數據B0curr按照如下方法計算：其中Δφcurr為當前場的相位差數據，按照如下方法計算：Δφcurr＝φ(r,tp2)curr-φ(r,tp1)curr其中，φ(r,tp1)curr為當前場第一相位數據，由EchoAcurr經快速傅里葉變換后計算相位得到；φ(r,tp2)curr為當前場第二相位數據，由EchoBcurr經快速傅里葉變換后計算相位得到。在本實施例中：本方法采集得到的當前場Profile如圖4-2所示，信噪比較圖4-1梯度回波采得的當前場數據有了明顯的提高。(3)根據步驟(1)和步驟(2)獲得的場圖矩陣和當前場分布數據，采用Tikhonov規(guī)則化的奇異值分解法，得到各勻場線圈的電流改變量。具體步驟為：采用Tikhonov規(guī)則化的奇異值(singularvaluedecomposition)分解法解如下方程組：FieldMap*x＝B0curr其中x為勻場線圈電流的改變量，將勻場線圈電流值修改x后，重新設置到勻場線圈中。在本實施例中：x對應7組勻場線圈的電流改變量，具體數值如下：)(4)根據步驟(3)中獲得的各勻場線圈的電流改變量，計算磁場的相位殘差：當相位殘差滿足收斂條件時結束，否則將各勻場線圈的電流修改后的磁場作為當前場，重復步驟(2)至(4)。所述收斂條件優(yōu)選為：相位殘差的L2范數小于預設的L2范數閾值，即||B0curr-FieldMap*x||2<ε其中，|B0curr-FieldMap*x|為相位殘差，ε為預設的L2范數閾值。在本實施例中：||B0curr-FieldMap*x||2＝0.18，ε＝0.2，梯度勻場收斂，結束梯度勻場。本方法與傳統梯度回波勻場方法的效果對比：如圖5所示，初始磁場條件下水樣在4.7ppm處1H譜的線性為97.0/725.3/1085.2(圖5-1)，相同的初始磁場均勻性、脈沖寬度、脈沖功率、采樣譜寬等實驗條件下，采用傳統的梯度回波自動勻場后采集單脈沖得到譜圖的線形為69.0/561.2/579.7(圖5-2)，而采用本方法梯度勻場后采集單脈沖得到譜圖的線形為51.5/191.7/203.8(圖5-3)。變換多組初始磁場條件，用傳統的梯度回波勻場方法和本方法進行對比實驗，結果如下：注：表中數據為單脈沖實驗得到譜圖的線性，分別表示譜峰在50％/0.55％/0.11％處的譜峰寬度。