基于區(qū)域內(nèi)容感知核范數(shù)的低劑量CT影像去噪

宋 蕓，張元科，盧虹冰，邢宇翔，馬建華

（1. 曲阜師范大學信息科學與工程學院，山東日照276800；2. 南方醫(yī)科大學生物醫(yī)學工程學院，廣州510515；3. 空軍軍醫(yī)大學生物醫(yī)學工程系，西安710032；4. 清華大學工程物理系，北京100084）

（?通信作者電子郵箱lemonzyk@fmmu.edu.cn）

0 引言

計算機斷層掃描（Computed Tomography，CT）產(chǎn)生的X 射線輻射暴露能夠誘發(fā)白血病、癌癥以及遺傳性等疾病，已在全球范圍內(nèi)引起極大關注［1-3］。為有效降低輻射劑量，臨床上最常采用操作是降低數(shù)據(jù)采集中X 射線管電流量（嚴格指管電流強度乘以曝光時間，單位為mAs），即低劑量CT（Low-Dose CT，LDCT）掃描。但LDCT 掃描引起的“光子饑餓”現(xiàn)象使得重建后的CT 影像存在嚴重的偽影噪聲［4］，影響診斷結(jié)果。另一方面，隨著先進影像分析技術的發(fā)展，“影像組學”（Radiomics）通過提取醫(yī)學影像中蘊含的大量高維特征，在臨床診斷與疾病預后方面取得突破性的進展［5］。而在影像組學中，影像的紋理細節(jié)特征對于分析至關重要。因此，迫切需要一種兼具噪聲抑制和紋理細節(jié)保持的LDCT影像去噪方法。

針對LDCT 影像中偽影噪聲，近年來學者們已提出大量去噪算法，按數(shù)據(jù)處理不同階段大致可分為三類:投影域濾波策略、迭代重建策略以及影像域濾波策略。投影域濾波策略直接對原始數(shù)據(jù)進行平滑濾波，并對濾波后數(shù)據(jù)采用傳統(tǒng)濾波反投影方法（Filtered Back Projection，F(xiàn)BP）方法重建；迭代重建策略［6］通過目標方程對原始數(shù)據(jù)及重建圖像先驗進行建模，并采用迭代優(yōu)化實現(xiàn)高質(zhì)量成像目的［6-7］；影像域濾波策略直接以FBP 重建影像為處理對象，根據(jù)圖像局部或非局部結(jié)構(gòu)特征，采用線性或非線性［8-15］方法進行去噪處理。

考慮到當前大部分商用CT 掃描設備僅提供重建后影像數(shù)據(jù)，為解決現(xiàn)實臨床應用問題，本文采用影像域濾波策略，目標是從含有噪聲的LDCT 影像y 估計其復原影像x，通?？杀硎緸閥 = x + e，其中e 表示加性噪聲。針對圖像中噪聲問題，當前去噪算法大致可以分為空間域去噪方法和變換域去噪方法?？臻g域方法主要是對影像像素灰度值進行處理，如非局部均值濾波、基于塊的局部最優(yōu)維納濾波［16-17］、全變分（Total Variation，TV）最小化方法［18-20］等；變換域方法是將影像從空間域轉(zhuǎn)換到變換域，通過對變換系數(shù)進行規(guī)整約束以達到去噪目的，如小波變換［21-22］、三維塊匹配濾波（Block Matching 3D，BM3D）［23］、主成分分析（Principle Component Analysis，PCA）變換［24-25］、字典學習［26-28］和變換學習［29］等。

近年來，利用影像的非局部自相似性進行去噪的方法取得了非常好的效果，成為該領域研究熱點。該方法通過相似影像塊聚類，利用每聚類組內(nèi)影像塊的相關性進行去噪。其中，基于低秩（Low Rank，LR）建模的算法表現(xiàn)出優(yōu)越的性能［30］。該方法主要包括四步:1）對于每個目標影像塊（尺寸為以列向量形式表示）yi∈Rn×1，利用塊匹配操作從y 中搜索與其相似的塊其中構(gòu)造數(shù)據(jù)矩陣其列是的列向量形式；3）計算Yi的低秩近似4）聚合（aggregate）中的影像塊（即影像塊放回其所在位置并對每像素上的所有影像塊取平均值），形成最終影像估計。其中，步驟3）通過最小化目標函數(shù)以獲取Yi的低秩近似解。這里R(Xi)是秩正則項，當前最常用方法是核范數(shù)最小化（Nuclear Norm Minimization，NNM）［31］。NNM 通過對矩陣Yi的每個奇異值進行同等強度約束，在此基礎上構(gòu)建易于計算的凸范數(shù)正則項。然而在實際求解中，由于每個奇異值均具有明確的物理意義，即較大奇異值通常表示圖像邊界輪廓等高能量特征，較小奇異值表示紋理細節(jié)等低能量特征，而其余極小奇異值則主要由噪聲產(chǎn)生。若對各奇異值進行同等強度約束，則忽略了圖像中各結(jié)構(gòu)成分的差異性，難以實現(xiàn)具有紋理保持的圖像去噪。此外，不同于自然影像，LDCT影像受原始數(shù)據(jù)采集及非線性FBP重建過程影響，其噪聲為非平穩(wěn)噪聲。傳統(tǒng)去噪算法通常并未考慮這一特性，而對影像各區(qū)域噪聲進行同強度處理，從而進一步加劇影像局部紋理丟失。

針對上述問題，本文提出一種具有區(qū)域內(nèi)容感知能力的加權核范數(shù)最小化LDCT 影像去噪算法——rwNNM（regionalcontent-aware weighted NNM）算法。該算法在NNM 基礎上，根據(jù)影像局部噪聲強度以及不同奇異值水平自適應設置核范數(shù)權重，從而有效提高核范數(shù)的局部信息表達能力，進而提高去噪后影像紋理細節(jié)保持能力。其中，為估計LDCT 影像局部噪聲，本文結(jié)合核范數(shù)優(yōu)化方法，通過利用矩陣奇異值分解（Singular Value Decomposition，SVD）與PCA 的關系，改進了傳統(tǒng)噪聲估計計算方法，極大提高了計算效率。此外，針對上述步驟1），本文在LDCT 影像局部噪聲估計基礎上，采用一種新的基于局部統(tǒng)計特性的相似影像塊搜索策略，進一步提高了影像塊匹配操作的準確性。針對臨床肺CT 成像的低劑量仿真數(shù)據(jù)及真實低劑量數(shù)據(jù)的實驗結(jié)果表明，該算法能夠在去除LDCT影像噪聲的同時，有效保持局部紋理細節(jié)信息。

1 本文算法rwNNM

其中:Yi= UΣVT是影像Yi的奇異值分解；Sβ/2(Σ)是在參數(shù)β/2下，對角矩陣Σ 的軟閾值函數(shù)，即對于Σ 中的每個對角元素Σjj，都有

傳統(tǒng)的NNM 中，所有奇異值被式（3）以閾值β/2進行同等強度收縮，忽略了LDCT 影像數(shù)據(jù)矩陣中各奇異值的差異性。而實際中不同奇異值代表影像中不同的結(jié)構(gòu)成分（如影像邊界輪廓、紋理細節(jié)及噪聲等），因此在設計去噪模型時，應根據(jù)不同奇異值進行自適應衰減處理，以實現(xiàn)有效去噪的同時保持紋理細節(jié)。顯然，傳統(tǒng)NNM 模型及其相應的求解方法不夠靈活，無法處理這些問題。

為了提高核范數(shù)的結(jié)構(gòu)表達能力，本文采用如下基于區(qū)域內(nèi)容感知的加權核范數(shù)最小化LDCT影像去噪模型:

其中:

w =[wi，1，wi，2，…，wi，N]T，wi，j＞0 為賦給λi，j(Xi)的非負權重。該模型通過設置權重向量以增強原始核范數(shù)的結(jié)構(gòu)表示能力。其中如何設置權重wi，j及設計相應求解算法是本文重點解決的問題。

1）權重設置。

在LDCT 影像中，較大奇異值表示人體組織結(jié)構(gòu)輪廓信息，較小奇異值表示各組織內(nèi)容紋理細節(jié)信息，而最后極小奇異值則可認為是主要由噪聲產(chǎn)生。因此在去噪過程中，應盡可能地保留較大的奇異值（即組織結(jié)構(gòu)信息）；同時需對較小的奇異值進行合理衰減（即去除噪聲同時需保持紋理細節(jié)）。為此，本文采用一種重加權機制來自適應地調(diào)整權重:

2）求解算法。

考慮到式（4）為非凸問題，利用封閉形式解決NNP 問題的子梯度方法［32］已不再適用。文獻［34］證明了可以通過加權奇異值軟閾值操作得到式（4）的閉合形式最優(yōu)解:

其中:Yi= UΣVT是數(shù)據(jù)矩陣Yi的奇異值分解；Sw/2(Σ)是在權重向量w 下對角矩陣Σ 的加權軟閾值函數(shù)，對于對角矩陣Σ中的每個對角元素Σjj，都有如下加權軟閾值函數(shù):

綜上所述，本文所提rwNNM算法流程如下。

下面將詳細解釋算法中所涉及的局部噪聲估計以及塊匹配操作。

3）局部噪聲估計。

由于LDCT 影像噪聲為非平穩(wěn)噪聲，因此對局部噪聲方差估計能夠有效提高算法區(qū)域結(jié)構(gòu)及紋理保持能力。為此本文采用文獻［24］中提出的基于PCA 的局部噪聲方差估計方法估計文獻［24］中，噪聲估計的思路是首先計算含噪目標影像塊矩陣的協(xié)方差矩陣并對矩陣進行SVD，然后利用SVD 系數(shù)進行噪聲近似。該方法需計算協(xié)方差矩陣，并對其進行SVD，這極大了增加了算法時間復雜性。為此，本文結(jié)合核范數(shù)優(yōu)化方法，利用數(shù)據(jù)矩陣SVD 及PCA 變換的關系，改進了噪聲估計計算方法（即利用數(shù)據(jù)矩陣Yi的SVD 直接估計，有效提高了計算效率。具體方法如下:

a）假設yi當前目標影像塊，通過塊匹配操作找到其非局部相似影像塊并組成矩陣Yi。

b）對Yi進行SVD:Yi= UΣVT。根據(jù)SVD 與PCA 關系（即Yi協(xié)方差矩陣公式Cˉ= YYT= UΣ2UT）可知，U 實際可看作Yi做PCA分解的基向量，而其分解系數(shù)即為αi= UTYi。

其中:γ為常量；median(·)表示取數(shù)據(jù)集的中值。

4）塊匹配方法。

傳統(tǒng)的塊匹配方法利用目標影像塊yi與待判別影像塊y′i間的歐氏距離d2(yi，y′i)，組成相似影像塊矩陣Yi。然而由于LDCT影像噪聲是非平穩(wěn)噪聲，傳統(tǒng)方法難以捕捉非局部影像塊之間的相關性，影響塊匹配精度。為解決這一問題，本文采用了統(tǒng)計近鄰（Statistical Nearest Neighbors，SNN）方法［35］實現(xiàn)塊匹配，該方法在歐氏距離的計算過程中，結(jié)合影像局部噪聲選取相似影像塊y′i的判別準則，即:

其中:ο是偏移參數(shù)。該方法詳細過程請參考文獻［35］。

2 實驗與結(jié)果分析

為驗證本文所提rwNNM 算法的有效性，對臨床肺CT 成像的低劑量仿真數(shù)據(jù)及真實低劑量臨床影像數(shù)據(jù)進行去噪處理，并與當前具有代表性的NNM 算法［31］、TV 算法［18］、變換學習（Transform Learning，TL）算法［29］進行比較。臨床患者均簽署書面知情同意書并通過倫理委員會審查。其中:臨床患者1 采用Siemens Sensation 16 CT 進行了常規(guī)劑量CT（Fulldose CT，F(xiàn)DCT）掃描（管電壓120 kVp，管電流100 mAs），該數(shù)據(jù)用于算法的定量仿真實驗；臨床患者2 采用Siemens Sensation 16 CT 分別進行了低劑量（管電壓120 kVp，管電流20 mAs）和常規(guī)劑量（管電壓120 kVp，管電流100 mAs）掃描（兩次掃描間隔1 個月），該患者數(shù)據(jù)用于算法的真實臨床數(shù)據(jù)實驗。兩位患者的其他掃描參數(shù)如下:機架旋轉(zhuǎn)一圈0.5 s，螺距1，校準16 × 0.75 mm，掃描層厚5 mm，重建層厚2 mm，重建層間隔1 mm，自動劑量控制功能關閉。文中所有影響均采用兩種臨床常用顯示窗:軟組織窗［864，1 264］HU和肺窗［0，1 000］HU。

2.1 參數(shù)設置

本文算法參數(shù)包括:δ 為迭代正則化參數(shù)，實驗中設為固定值0.1；γ為噪聲估計參數(shù)，根據(jù)文獻［34］建議，實驗中設為固定值0.8；ο 為塊匹配參數(shù)，根據(jù)文獻［35］建議，實驗中設為固定值1；K 為迭代次數(shù)，實驗中設為10；經(jīng)反復實驗，本文算法影像塊尺寸為8× 8，塊匹配中非局部相似影像塊數(shù)量為150。其他算法均根據(jù)其相關文獻建議選取能夠取得最佳效果的參數(shù)。

2.2 計算機仿真實驗

2.2.1 視覺效果

實驗中采用患者1 的常規(guī)劑量CT 影像數(shù)據(jù)（如圖1（a））進行數(shù)字仿真:首先采用扇形投影并采用上述掃描參數(shù)對常規(guī)劑量影像進行正投影，然后根據(jù)文獻［36］方法在投影數(shù)據(jù)中添加噪聲（仿真劑量約為20 mAs），最后進行FBP 重建，如圖1（b）所示，從中可發(fā)現(xiàn)LDCT 影像中存在明顯的條帶狀非平穩(wěn)偽影噪聲。對比圖1 中不同算法的處理結(jié)果可發(fā)現(xiàn):NNM、TV 和TL 算法盡管能抑制部分噪聲，但是在軟組織區(qū)域和肺區(qū)域仍存在較強的條狀偽影噪聲，且影像紋理保持較差；本文rwNNM 算法能夠在有效去噪的同時保存紋理細節(jié)信息，尤其是矩形實線框線內(nèi)含肺結(jié)節(jié)感興趣區(qū)。

圖1 仿真實驗結(jié)果Fig.1 Simulation results

2.2.2 定量分析

本文采用法向量流（Normal Vector Flow，NVF）［37］、結(jié)構(gòu)相似指標（Structural SIMilarity index，SSIM）［38］和均方根誤差（Root Mean Square Error，RMSE）對所提rwNNM 算法進行定量分析。其中，NVF 指標評估各算法處理結(jié)果與參考標準FDCT 影像之間的紋理相似性。若算法NVF 中的箭頭有序并且與參考影像的NVF 方向一致，表明該算法的細節(jié)紋理保存性能較好。圖2 為各算法在肺癌結(jié)節(jié)區(qū)域（圖1（a）中虛線框內(nèi)區(qū)域）的NVF 影像，從中可看出:相對于NNM、TV 和TL 算法，本文rwNNM 算法所產(chǎn)生的NVF 箭頭與參考標準影像的NVF 匹配程度最高，說明rwNNM 算法能夠更好地保持影像的細節(jié)紋理特征。

本文進一步引入SSIM 和RMSE 對所提rwNNM 算法進行定量分析。其中，RMSE 是表征實驗結(jié)果準確性的度量指標，而SSIM 則是表征影像結(jié)構(gòu)保持性能的度量指標。所分析區(qū)域為包含肺結(jié)節(jié)的局部感興趣區(qū)。SSIM和RMSE的計算公式如下:

其中:r 和x 分別代表參考仿真模型和實驗結(jié)果的ROI I 區(qū)域；rˉ和xˉ是每個ROI I 區(qū)域的平均強度；σr和σx是標準偏差；σxr為參考仿真模型和實驗結(jié)果之間的協(xié)方差；I 代表ROI I 區(qū)域中像素的個數(shù)；c1=（K1Ls）2和c2=（K2Ls）2為常量，Ls為影像動態(tài)范圍。根據(jù)文獻［38］，K1和K2分別設置為0.01 和0.03。各算法的SSIM 和RMSE 定量分析結(jié)果如圖3 所示。從圖3 可看出:本文rwNNM 算法具有最低的RMSE 指標和最高的SSIM指標（兩指標分別提高14.38%和8.75%以上），表明本文算法在噪聲抑制及結(jié)構(gòu)相似性保持定量指標上優(yōu)于其他算法。

圖2 肺結(jié)節(jié)區(qū)域的NVF圖Fig.2 NVF images of lung nodule region

圖3 不同算法SSIM和RMSE分析結(jié)果Fig.3 Analysis results of RMSE and SSIM of different algorithms

2.3 臨床數(shù)據(jù)實驗

本節(jié)采用患者2 的120 kVp/20 mAs 低劑量CT 影像進行算法臨床LDCT 數(shù)據(jù)驗證。并用該患者一個月后常規(guī)劑量掃描數(shù)據(jù)作為結(jié)果評判參考。所有影像均采用兩種顯示窗:軟組織窗［864，1 264］HU和肺窗［0，1 000］HU。

從如圖4 所示的不同算法處理結(jié)果中可看出:本文rwNNM 算法在處理真實低劑量CT 數(shù)據(jù)時較其他算法仍能取得更好的噪聲抑制以及紋理細節(jié)保持效果。

為進一步對rwNNM 算法的臨床效果進行定性評估，本文取患者2 的15 個斷層采用各算法進行處理，并邀請三位具有5 年以上臨床經(jīng)驗的放射醫(yī)生對不同算法得到的影像結(jié)果（從降噪/分辨率/紋理保持三個方面）進行評分，分數(shù)范圍從0到10（0為最差，10為最好）。為保證客觀性，評分過程通過在電腦上隨機播放不同算法得到的去噪圖像，由每位醫(yī)生分別評分。最后計算每個算法的平均得分及標準差，如表1 所示。表1結(jié)果表明，從臨床放射醫(yī)生角度，rwNNM 算法較NNM、TV和TL算法能夠取得更好的臨床影像處理效果，且具有顯著統(tǒng)計差異性（P＜0.05）。

表1 放射醫(yī)生對不同算法處理結(jié)果的評分（均值±標準差）Tab. 1 Radiologists'scores of the results obtained by different algorithms（mean value±standard deviation）

圖4 臨床實驗結(jié)果Fig.4 Clinical results

3 結(jié)語

本文提出了具有區(qū)域內(nèi)容感知能力的加權核范數(shù)最小化LDCT 影像去噪算法，該算法針對低劑量CT 影像的非平穩(wěn)噪聲特點，采用基于局部統(tǒng)計特性的相似影像塊搜索策略，根據(jù)影像局部噪聲強度以及不同奇異值水平自適應設置核范數(shù)權重，有效提高核范數(shù)的局部信息表達能力，從而提高LDCT 復原后影像的紋理細節(jié)保持能力，進而有效提高基于LDCT 影像的臨床診斷及影像學分析的準確性。但是該算法在高噪聲強度下，LDCT 影像紋理細節(jié)保持能力下降。在未來的研究中，將進一步結(jié)合LDCT 投影數(shù)據(jù)的統(tǒng)計特性以及影像局部特性和非局部特性，開展基于區(qū)域內(nèi)容感知能力的加權核范數(shù)最小化與變換學習的LDCT影像重建相關研究。