本發(fā)明涉及一種核磁共振信號處理技術(shù)，特別涉及一種基于PDCO的低場核磁共振二維譜反演算法。

背景技術(shù)：

在世界范圍內(nèi)，核磁共振技術(shù)發(fā)展迅猛，在很多領(lǐng)域得到了廣泛的發(fā)展，比如高場強(qiáng)核磁共振技術(shù)應(yīng)用于對人體的臨床診斷等。隨著核磁共振技術(shù)的不斷發(fā)展，人們發(fā)現(xiàn)核磁共振技術(shù)不僅可以應(yīng)用于臨床，還可以在其他領(lǐng)域(如食品科學(xué)、農(nóng)業(yè)、石油能源、材料科學(xué)、紡織化工等)發(fā)揮其他科學(xué)儀器所不能發(fā)揮的作用。很多科學(xué)研究已經(jīng)證明，在這些領(lǐng)域中應(yīng)用核磁共振分析技術(shù)，可以解決現(xiàn)有的其他科學(xué)儀器所不能解決的問題，對這些領(lǐng)域的科學(xué)進(jìn)步起到了舉足輕重的推動作用。與用于臨床的高場核磁共振技術(shù)不同的是，在食品科學(xué)、農(nóng)業(yè)、石油能源、材料科學(xué)、紡織化工等領(lǐng)域采用低場核磁共振技術(shù)就可以很好的解決相關(guān)問題。此外，低場核磁共振分析儀器由于采用低場強(qiáng)磁體，可以使所研發(fā)的分析儀器體積較小，安裝、調(diào)試、維護(hù)、操作都很方便。因此，低場核磁共振分析儀器得到了科學(xué)界的關(guān)注。

CPMG(Carr-Purcell-Meiboom-Gill)序列速度快，是低場核磁共振中最常用的序列之一。研究者常常利用CPMG等序列的原始數(shù)據(jù)和樣品橫向弛豫時間、縱向弛豫時間的分布特點，進(jìn)行時域譜反演的相關(guān)研究。但是核磁共振采集到的原始信號并不能得到直觀的樣品結(jié)構(gòu)信息，需要通過反演技術(shù)才能得到易于理解的時域譜。一般，一維時域譜就可以為研究者分析樣品的組成、性質(zhì)等提供重要的依據(jù)。但是，隨著低場核磁共振應(yīng)用的不斷深入，研究者們發(fā)現(xiàn)實驗得到的一維譜可能存在峰的交疊。解決這一問題通常需要進(jìn)行多次額外實驗或者引入其它輔助手段，二維時域譜技術(shù)應(yīng)運而生。二維譜不僅簡化了原有流程，還可以為分析樣品提供更多有價值的信息。

國內(nèi)二維反演技術(shù)由于起步晚，國外相關(guān)企業(yè)的技術(shù)壟斷，加上二維譜反演需要處理的數(shù)據(jù)量巨大等問題，想要發(fā)展具有極大的挑戰(zhàn)性。目前二維反演算法，國外研究者以顯式增加罰函數(shù)的正則化算法為主，對解的某些性質(zhì)進(jìn)行限制，以得到合理的反演結(jié)果；國內(nèi)研究者最早使用基于截斷奇異值分解的二維反演算法，并取得了很多成果，在信噪比(Signal to Noise Ratio,SNR)較高的情況下能夠得到可靠的反演結(jié)果。現(xiàn)在一般采用的是標(biāo)準(zhǔn)Tikhonov正則化形式，懲罰項用于限制解的范數(shù)，這種方式的問題求解策略也相對成熟，此外還有對斜率和曲率進(jìn)行平滑的正則化形式。這些正則化算法的關(guān)鍵在于如何選擇一個合適的正則化因子，現(xiàn)有正則化因子的選擇算法需要一定的人為干預(yù)，還待完善。

技術(shù)實現(xiàn)要素：

本發(fā)明是針對核磁共振信號處理中反演技術(shù)存在的問題，提出了一種低場核磁共振二維譜反演算法，在PDCO(a Primal-Dual interior method for Convex Objectives一種求解凸規(guī)劃的原始對偶內(nèi)點法)基礎(chǔ)上推廣到二維空間，用于求解二維反演問題。PDCO迭代使用L1范數(shù)與L2范數(shù)，使得反演結(jié)果更貼近于真實的譜分布，并且保證了解的非負(fù)性，不再需要額外的非負(fù)約束計算。該算法能夠提高反演運算的計算精度和魯棒性，可以得到清晰、準(zhǔn)確的二維譜。

本發(fā)明的技術(shù)方案為：一種低場核磁共振二維譜反演算法，具體包括如下步驟：

1)讀取低場核磁共振設(shè)備采集得到的原始數(shù)據(jù)文件，提取出數(shù)據(jù)文件中包含的采樣時間和對應(yīng)時刻的采樣數(shù)據(jù)M；

2)對原始數(shù)據(jù)進(jìn)行橫向弛豫時間、縱向弛豫時間布點，矩陣的拼接、張量積預(yù)處理操作得到反演核心矩陣K和信號幅值m；

3)計算采集數(shù)據(jù)的信噪比，根據(jù)信噪比計算反演參數(shù)，采用求解凸規(guī)劃的原始對偶內(nèi)點法PDCO解決反演問題表達(dá)式如下：

${\min }_{s,r}{\mathrm{\λ}}_1\left|\right|s|{|}_1+\frac{1}{2}{\mathrm{\λ}}_2|\left|s\right|{|}_2^2+\frac{1}{2}\left|\right|r|{|}_2^2$

s.t.Ks+r＝m

s≥0

其中s＝vect(S)，vect表示將矩陣按列拼接，形成一個列向量，S(T₂，T₁)表示橫向弛豫時間為T₂、縱向弛豫時間為T₁的物質(zhì)的含量，反演參數(shù)λ₁＝α₁||m||₁/SNR,反演參數(shù)λ₂＝α₂/SNR，r表示擬合殘差，||*||₁表示取某向量的L1范數(shù)，||*||₂表示取某向量的L2范數(shù)，信噪比SNR使用的計算方法為：采樣數(shù)據(jù)信號最大值除以衰減后數(shù)據(jù)的方差，α₁與α₂參數(shù)值由大量仿真實驗測得：α₁取值為0.0001，α₂為20；

4)進(jìn)行迭代計算，每次迭代并保證解的非負(fù)性，得到最優(yōu)解s'；

5)生成網(wǎng)格，繪出二維譜，用T₂、T₁這兩組數(shù)據(jù)分別表示Y軸X軸信息，將這兩組信息生成網(wǎng)格，根據(jù)T₂、T₁中元素的個數(shù)并對最優(yōu)解s'進(jìn)行重新排序成s”，s”表示二維譜反演結(jié)果信號的幅值，利用T₂、T₁、s”這3組數(shù)據(jù)繪出二維譜。

本發(fā)明的有益效果在于：本發(fā)明一種低場核磁共振二維譜反演算法，與目前國內(nèi)外文獻(xiàn)報道的其他低場核磁共振反演算法相比，其使用L1正則化和L2正則化而不是通常的只使用L2正則化，更貼近原始信號的稀疏性，得到的解更能體現(xiàn)真實的譜分布；原始對偶內(nèi)點法在計算時就保證了解的非負(fù)性，不需要進(jìn)行額外的非負(fù)約束計算，簡化了反演流程；對緊密相連的臨近峰也能夠區(qū)分，計算精度高；魯棒性好，在不同信噪比數(shù)據(jù)情況下，能夠得到穩(wěn)定的反演結(jié)果。

附圖說明

圖1為本發(fā)明低場核磁共振二維譜反演算法步驟示意圖；

圖2-1為本發(fā)明信噪比為1000的構(gòu)造高斯峰二維譜仿真實驗結(jié)果圖；

圖2-2為本發(fā)明信噪比為1000的CPMG數(shù)據(jù)串仿真實驗結(jié)果圖；

圖2-3為本發(fā)明信噪比為1000的反演仿真實驗結(jié)果圖；

圖2-4為本發(fā)明信噪比為1000的標(biāo)準(zhǔn)Tikhonov正則化算法的反演仿真實驗結(jié)果圖；

圖3-1為本發(fā)明信噪比為100的構(gòu)造高斯峰二維譜仿真實驗結(jié)果圖；

圖3-2為本發(fā)明信噪比為100的CPMG數(shù)據(jù)串仿真實驗結(jié)果圖；

圖3-3為本發(fā)明信噪比為100的反演仿真實驗結(jié)果圖；

圖3-4為本發(fā)明信噪比為100的標(biāo)準(zhǔn)Tikhonov正則化算法的反演仿真實驗結(jié)果圖；

圖4-1為本發(fā)明實驗案例低濃度溶液反演結(jié)果圖；

圖4-2為本發(fā)明實驗案例加入高濃度溶液反演結(jié)果圖。

具體實施方式

下面結(jié)合附圖對本發(fā)明提供的具體實施方式作詳細(xì)說明。

如附圖1所示，一種基于PDCO的低場核磁共振二維譜反演算法，包括如下步驟：

a.讀取低場核磁共振設(shè)備采集得到的原始數(shù)據(jù)文件。

讀取低場核磁共振設(shè)備采集得到的原始數(shù)據(jù)文件，提取出數(shù)據(jù)文件中包含的采樣時間和對應(yīng)時刻的采樣數(shù)據(jù)M等信息。

b.對原始數(shù)據(jù)進(jìn)行橫向弛豫時間、縱向弛豫時間布點，矩陣的拼接、張量積等預(yù)處理操作得到反演核心矩陣K和信號幅值m。

二維反演問題就是求解如式(1)所示的具有兩個核的Fredholm積分方程，τ₂表示各回波波峰的回波時刻，τ₁表示等待時間，S(T₂，T₁)表示橫向弛豫時間為T₂、縱向弛豫時間為T₁的物質(zhì)的含量，M表示特定時刻的采樣數(shù)據(jù)。

$M({\mathrm{\τ}}_2,{\mathrm{\τ}}_1)={\∫}_{T_{2\min }}^{T_{2\max }}{\∫}_{T_{1\min }}^{T_{1\max }}{e}^{-\frac{{\mathrm{\τ}}_2}{T_2}}\·S(T_2,T_1)\·(1-2e^{-\frac{{\mathrm{\τ}}_1}{T_1}}){\mathrm{dT}}_1{\mathrm{dT}}_2---\left(1\right)$

將公式(1)中的弛豫時間、離散化后可以用矩陣的形式表示如式(2)所示，其中(上標(biāo)i,j表示矩陣中元素位置的索引，TE表示CPMG序列中的回波時間，TW為等待時間)。

$M=K_2{\mathrm{SK}}_1^T---\left(2\right)$

上標(biāo)T表示矩陣轉(zhuǎn)置，對于公式(2)，可以利用矩陣的張量積將兩個核K₁、K₂合成一個核，合并后公式如(3)式所示，其中m＝vect(M)，s＝vect(S)，vect表示將矩陣按列拼接，形成一個列向量，表示兩個矩陣的張量積。

m＝Ks (3)

這樣轉(zhuǎn)換之后，二維反演問題就變成一個已知m和K求s的問題，其中矩陣K稱為反演的核心矩陣。

為了得到符合真實情況的譜分布，要對T₂、T₁進(jìn)行布點，布點方法是在以10為底的對數(shù)空間內(nèi)均勻布點，布點數(shù)目根據(jù)具體實驗情況而定，實驗精度要求越高，布點數(shù)目越多，一般為了節(jié)約計算時間將布點數(shù)目設(shè)為64。

該步驟中所述的數(shù)據(jù)預(yù)處理主要包括：T₂、T₁布點；按照公式(2)中K₁、K₂的定義計算得到K₁、K₂；利用矩陣的張量積轉(zhuǎn)化K₁、K₂得到K；使用矩陣的拼接方法得到信號幅值m等操作。

c.計算采集數(shù)據(jù)的信噪比，根據(jù)信噪比計算反演參數(shù)。PDCO解決上述反演問題表達(dá)式如下：

${\min }_{s,r}{\mathrm{\λ}}_1\left|\right|s|{|}_1+\frac{1}{2}{\mathrm{\λ}}_2|\left|s\right|{|}_2^2+\frac{1}{2}\left|\right|r|{|}_2^2---\left(4\right)$

s.t.Ks+r＝m

s≥0

其中反演參數(shù)λ₁＝α₁||m||₁/SNR,反演參數(shù)λ₂＝α₂/SNR，r表示擬合殘差，*||₁表示取某向量的L1范數(shù)，||*||₂表示取某向量的L2范數(shù)。信噪比SNR使用的計算方法為：采樣數(shù)據(jù)信號最大值除以衰減后數(shù)據(jù)的方差。α₁與α₂參數(shù)值由大量仿真實驗測得：α₁取值為0.0001，α₂為20(信噪比較高時，可將λ₂直接設(shè)為0.5)。

d.進(jìn)行迭代計算，每次迭代并保證解的非負(fù)性，得到最優(yōu)解。

迭代計算時主要求解如下問題:

${\min }_{s,r,s_1}{\mathrm{\λ}}_1\left|\right|s|{|}_1+\frac{1}{2}{\mathrm{\λ}}_2|\left|s\right|{|}_2^2+\frac{1}{2}\left|\right|r|{|}_2^2-\mathrm{\μ}\underset{l=1}{\overset{n}{\Σ}}ln\left({\[s_1\]}_l\right)---\left(5\right)$

s.t.Ks+r＝m :y

s-s₁＝0 :z₁

其中μ,s₁≥0(μ表示懲罰因子，通常取一個很小的正參數(shù)，s₁表示松弛變量，通常第一次迭代時所有元素初始化為1)，l表示s₁中元素的位置，n表示s₁中元素總個數(shù)，由此實現(xiàn)累加求和。引入懲罰因子與松弛變量，保證了PDCO算法能夠拓展到二維空間，并且增加了對數(shù)項保證了迭代過程中解的非負(fù)性。問題轉(zhuǎn)換之后，在求解時將上式與Karush-Kuhn-Tuucker(KKT)條件聯(lián)立，使用牛頓法進(jìn)行搜索，即可進(jìn)行迭代計算。閾值誤差設(shè)置為0.001，達(dá)到終止條件時，停止迭代，得到最優(yōu)解s'。

e.生成網(wǎng)格，繪出二維譜。用T₂、T₁這兩組數(shù)據(jù)分別表示Y軸X軸信息，將這兩組信息生成網(wǎng)格。根據(jù)T₂、T₁中元素的個數(shù)并對最優(yōu)解s'進(jìn)行重新排序成s”，s”表示二維譜反演結(jié)果信號的幅值，利用T₂、T₁、s”這3組數(shù)據(jù)繪出二維譜。

本發(fā)明的效果能夠通過以下實驗進(jìn)一步說明。

1.仿真實驗：

實驗首先構(gòu)造一個中心位于(T₁，T₂)＝(100,10)ms處的高斯峰作為理想的T1-T2譜，然后向正演結(jié)果中加入一定程度的高斯白噪聲，得到不同信噪比的仿真數(shù)據(jù)。

2.仿真實驗結(jié)果及結(jié)果分析

圖2和圖3分別是信噪比為1000(仿真低場下的高信噪比數(shù)據(jù))和信噪比為100(仿真低場下的低信噪比數(shù)據(jù))的仿真結(jié)果，這兩幅圖中2-1/3-1圖顯示的是構(gòu)造的高斯峰二維譜；2-2/3-2圖是原始數(shù)據(jù)正演得到的8條CPMG數(shù)據(jù)串；2-3/3-3圖是本發(fā)明提出的基于PDCO的二維譜反演算法反演結(jié)果；2-4/3-4圖是使用標(biāo)準(zhǔn)Tikhonov正則化算法的反演結(jié)果。

將這兩組實驗進(jìn)行對比分析，可以看出在高信噪比環(huán)境下，2種算法都能得到反演結(jié)果，這兩種算法都有一定的干擾峰出現(xiàn)，但明顯可以看出基于PDCO的二維反演算法結(jié)果優(yōu)于傳統(tǒng)的標(biāo)準(zhǔn)Tikhonov正則化算法。在低信噪比環(huán)境下，傳統(tǒng)的Tikhonov正則化算法，反演結(jié)果比較穩(wěn)定，譜峰仍然明顯被拓寬，有虛假輪廓出現(xiàn)；基于PDCO的混合反演算法得到的反演結(jié)果中也產(chǎn)生了虛假輪廓，反演結(jié)果相對于另一種算法明顯更貼近仿真數(shù)據(jù)。

為了進(jìn)一步驗證本發(fā)明算法的抗噪性能，表1列舉了信噪比分別為3000、1000、300、200、100、50下的兩種算法對擬合誤差比較結(jié)果。由表1可知，本發(fā)明的算法與標(biāo)準(zhǔn)Tikhonov正則化算法能夠在不同信噪比數(shù)據(jù)的情況下得到穩(wěn)定的結(jié)果；本發(fā)明的算法與標(biāo)準(zhǔn)Tikhonov正則化算法相比，擬合誤差更小，反演結(jié)果更貼近真實的仿真數(shù)據(jù)。

表1

3.實驗案例：

本實驗使用NIUMAG公司的NMI-20低場核磁共振分析儀進(jìn)行采樣，采樣序列為IR-CPMG，TW設(shè)置為0～2500ms內(nèi)對數(shù)均勻分布的8個點(對應(yīng)界面參數(shù)DL1)，TE設(shè)為0.235ms，累加次數(shù)NS為4。實驗樣品為預(yù)先制備的兩種不同濃度的CuSO4溶液，分別密封在兩個無核磁信號的色譜瓶中。

4.實驗案例結(jié)果及結(jié)果分析：

實驗首先在無磁試管中放入裝有低濃度溶液的色譜瓶進(jìn)行采樣，使用本發(fā)明算法反演得到了圖4-1所示的T1-T2譜；然后再向試管中放入裝有高濃度溶液的色譜瓶(色譜瓶足夠小，保證樣品均處于設(shè)備的勻場區(qū)域)，得到了圖4-2所示的T1-T2譜。

通常，在水中添加無磁性重金屬離子會影響質(zhì)子的弛豫時間，添加的離子濃度越高，采集到的CPMG信號衰減越快，質(zhì)子的弛豫時間越短，在二維譜中的位置就越靠近原點。該實驗中，加入高濃度的CuSO4溶液后，在T1-T2譜中更靠近原點的位置又出現(xiàn)了一個譜峰，說明樣品中出現(xiàn)了更短弛豫時間的成分，進(jìn)而證明了無磁性重金屬離子對質(zhì)子影響的論斷。