潘玲佼，范偉偉，吳全玉

（1.江蘇理工學院電氣信息工程學院生物信息與醫(yī)藥工程研究所，常州 213001；2.常州大學微電子與控制工程學院，常州 213164）

引 言

光學相干層析成像（Optical coherence tomography，OCT）技術是1991年問世的一種基于低相干干涉的具有高縱向空間分辨率的光學成像技術。它具有非接觸性、無創(chuàng)生物組織成像等優(yōu)勢，其最重要的應用之一就是眼科疾病的篩查和診療［1］。隨著2007年商業(yè)化譜域OCT（Spectral?domain OCT，SD?OCT）技術的問世，實現(xiàn)了三維高分辨率成像，從而進一步推動了OCT 技術在眼科的應用，為青光眼、年齡相關性黃斑變性和糖尿病黃斑水腫等視網膜疾病的定量分析和研究提供了重要幫助［2?6］。在SD?OCT 中，當OCT 系統(tǒng)對樣品組織探測時，將樣品組織等效為一系列的反射面，則可探測某深度上組織的反射率，形成A 掃描成像。將光束對組織樣品進行橫向移動，得出多個A 掃描圖像，拼接出一幅二維斷層圖像，即B 掃描圖像。如圖1 所示，三維SD?OCT 成像由若干二維B 掃描圖像組成，每個B 掃描圖像又由若干A 掃描組成。OCT 成像系統(tǒng)采用了干涉技術，受低相干光源和成像樣本自身結構的共同影響，OCT 圖像中含有大量散斑。散斑的相關特性可用來描述微血管系統(tǒng)，反映血流速度等信息，對診斷有一定的價值。然而，由于散斑是多個散射元疊加的結果，在B 掃描圖像中表現(xiàn)為大量隨機高反差亮斑，這些亮斑稱為散斑噪聲。散斑噪聲淹沒了樣本的結構細節(jié)，導致圖像質量退化，嚴重影響人們對圖像信息的準確獲取，給醫(yī)生進行疾病診斷以及計算機輔助疾病診斷分析帶來了極大挑戰(zhàn)。如何抑制散斑噪聲一直是國內外學者致力研究的重要課題。

圖1 三維光學相干層析成像示意圖Fig.1 Illustration of 3D optical coherence tomography

目前，針對散斑去噪主要有硬件和軟件兩種處理方法［7?9］。硬件方法主要依靠改進光學系統(tǒng)架構，該方法能夠有效抑制散斑噪聲，但也會使得OCT 成像系統(tǒng)結構更加復雜，造價更高。與硬件方法相比，軟件方法無需改變系統(tǒng)結構，更加靈活、易實現(xiàn)。商用OCT 掃描儀提升OCT 成像質量的通用方法是在樣本同一位置重復掃描多次，將多次掃描圖像疊加取平均來消除散斑噪聲。重復掃描圖像數(shù)量越多，則獲得的B 掃描的成像質量越高。然而，隨著重復掃描圖像數(shù)量的增加，相應的掃描時間也會延長?？紤]到長時間掃描過程中，眨眼等微小動作都會嚴重影響成像的質量，因此，需要限制掃描時間。在掃描時間有限的前提下，能完成的掃描次數(shù)也是有限的。對同一位置多次重復掃描能有效提高單張B 掃描的成像質量，但會相應減少三維OCT 中B 掃描的數(shù)量，從而降低三維成像分辨率。而增加B 掃描的數(shù)量又勢必降低單張B 掃描的成像質量。本文擬用低秩矩陣理論來解決這一問題。現(xiàn)有的抑制散斑噪聲的方法，如空間域局部統(tǒng)計濾波算法［10］、各項異性擴散濾波算法［11］、基于小波變換的濾波算法［12］、基于塊匹配的方法［13］等，盡管能在一定程度上有效去除散斑噪聲，但并不適合應用在商用OCT 掃描儀上。近年來也有學者嘗試用深度學習方法解決散斑問題［14］，然而該方法需要大量圖像樣本進行訓練，且訓練時間較長。 王帆等［15］提出了使用主成分分析（Principle component analysis，PCA）方法估計OCT 圖像散斑噪聲水平。它將OCT 圖像原始信號矩陣進行奇異值分解，然后利用一個或幾個較大奇異值估計去噪圖像。與其他方法相比，該方法僅利用少量相鄰幀生物結構的高度相似性，能夠在確保去除散斑噪聲的同時較好地保留原有生物結構信息。但是這種方法如果對主奇異值選取不當，則會造成散斑噪聲估計不準確且魯棒性較差。袁治靈等［16］提出了一種基于穩(wěn)健性主成分分析算法（Robust PCA，RPCA）的光學相干層析成像散斑噪聲去除方法，將OCT 原始圖像分解為散斑噪聲圖像和樣品截面圖像的最佳估計之和，一定程度上解決了PCA 方法的缺點。但是該方法只能應用于具有唯一低秩稀疏分解結構的矩陣，并且對于高維數(shù)據(jù)，求解運算復雜度較高，對于噪聲的適應性較差、分解效率低。本文在前人研究基礎上，考慮視網膜OCT 成像中眼球運動帶來的運動偽差特性，結合雙邊隨機投影算法，提出了一種新的光學層析成像去噪方法。該方法使用雙邊隨機投影代替穩(wěn)健性主成分分析，降低了計算的時間復雜度，提高了分解速度。通過優(yōu)化投影可以精確恢復復雜信號，從而準確提取目標信號。

1 基本方法

1.1 算法概述

在三維視網膜OCT 掃描圖像中，由于連續(xù)掃描的相鄰二維圖像的生物組織結構信息幾乎相同，利用這種組織結構之間高度相似性以及OCT 圖像高分辨率特性，可將三維OCT 掃描圖像中相鄰的若干二維斷層圖像（即B 掃描圖像）中相似的生物組織結構看作一個低秩矩陣。因此，OCT 散斑去噪問題即等同于從被污染的信號中恢復原始低秩矩陣的問題。在視網膜OCT 圖像掃描中，信號污染主要來自于兩方面，一方面是由于圖像掃描過程中人的呼吸、心跳等不可避免因素引起眼球抖動，從而造成了視網膜OCT 圖像的運動偽差［17］；另一方面主要來自于掃描中不可避免的散斑噪聲。運動偽差信號是少數(shù)異常體，相對比較稀疏，具有一定稀疏度，因此可看作是一個稀疏矩陣。而散斑噪聲信號則看作是一個噪聲矩陣。

圖2 給出了算法流程圖。首先對OCT 圖像進行對齊預處理，增強相鄰幀之間的相關性；其次將原始OCT 圖像信號分解為無噪低秩矩陣部分、對應運動偽差的稀疏矩陣部分以及噪聲矩陣3 部分；最后采用雙邊隨機投影算法求解，提取其中目標信號，從而達到去除噪聲、恢復無噪圖像的目的。

圖2 算法流程圖Fig.2 Flow chart of the pro?posed algorithm

1.2 OCT 圖像預處理

如圖1 所示，對每個OCT 圖像，在慢掃描方向，即Y方向（y?z視圖），相鄰兩個B 掃描圖像之間存在很大失真，主要表現(xiàn)為相鄰兩個B 掃描圖像跳變。這主要是由于眼睛和頭部在垂直方向運動引起的。因此在去噪之前，為保證相鄰幀生物組織結構之間的高度相似性，首先需進行對齊處理，從而保證信號提取的準確性。

B 掃描圖像對齊步驟如下。首先分割出OCT 圖像中視網膜上邊界，因為該邊界在視網膜所有邊界中最為突出，比較容易準確地進行分割，本文采用多分辨率表面檢測方法［18］進行分割。分割后計算B掃描圖像中視網膜上邊界的平均深度。假設M×N大小的第j幀B 掃描圖像分割后視網膜上邊界各像素點的深度分別為dij（i=0，1，2，…，N-1），則第j幀B 掃描圖像視網膜上邊界的平均深度為

由于視網膜上邊界中心凹是自然凹的曲面，因此，式（1）只使用除中心凹外左側和右側各t個像素點，t≈0.2N。以最小作為基準dbase，則各B 掃描圖像的位移值為

各B 掃描圖像根據(jù)Δdj縱向位移進行對齊，對齊后所有B 掃描圖像視網膜上邊界的平均深度一致。

視網膜上表面校正前后效果如圖3 所示。對齊之前相鄰B 掃描之間有很多跳變，視網膜上表面高低起伏。對齊之后相鄰B 掃描之間平滑連續(xù)，視網膜上表面也趨于平滑。對齊之后大大減少了運動偽差的影響，視網膜各層趨于平滑，視網膜成像形狀更為逼真。

圖3 OCT 圖像預處理前后的三維渲染圖Fig.3 3D rendering of the OCT image before and after preprocessing

1.3 OCT 信號分解

根據(jù)低秩分解理論，退化圖像矩陣X可分解為低秩矩陣L、稀疏矩陣S和噪聲矩陣N之和，退化前的數(shù)據(jù)可通過低秩矩陣來逼近。

在三維視網膜OCT 掃描圖像中，由于連續(xù)掃描的相鄰二維圖像的生物組織結構信息幾乎相同，若將一幀包含m個像素的二維掃描圖像看作矩陣列向量，相鄰n幀二維掃描圖像中相似的生物組織結構則可看作低秩矩陣，而該n幀圖像中的運動偽差和噪聲可分別看作稀疏矩陣和噪聲矩陣。OCT 設備采集的相鄰n幀數(shù)據(jù)則可看作退化圖像矩陣。本文取n為5 以保證生物組織結構的相似度，其矩陣關系式為

抑制OCT 圖像噪聲即要使N能量最小，即

式中：‖ ?‖F(xiàn)為Frobenius 范數(shù)，rank 為矩陣秩，card 為矩陣分量個數(shù)。由于式（5）中L和S均為未知矩陣，給定r和k，式（5）的優(yōu)化問題可轉化為以下兩個交替求解直至收斂的子問題，即

式中：Lt和St分別為t次迭代后的低秩部分和稀疏部分。

1.4 雙邊隨機投影散斑去噪

盡管奇異值分解的低秩逼近能夠達到對原始數(shù)據(jù)的最小重構誤差，但這一方法用于三維OCT 圖像這樣的高維數(shù)據(jù)計算復雜度較高。本文采用了雙邊隨機投影低秩逼近方法來近似和加速奇異值分解低秩逼近。根據(jù)雙邊隨機投影算法［19］，若有矩陣X∈Rm×n(m>n) 的雙邊隨機投影Y1=XA1和Y2=XTA2，A1∈Rm×r，A2∈Rn×r均為高斯隨機矩陣，則可構造X的r秩低秩矩陣L，即

為進一步提升低秩逼近的精度，使用power scheme 模型［20］進行優(yōu)化。用矩陣=(XXT)q X替代X計算雙邊隨機投影。因為與X相比，的奇異值衰減更快，且X與的奇異向量相同。的雙邊隨機投影為

為進一步優(yōu)化低秩逼近的速度和精度，對Y1和Y2進行矩陣正交三角分解（QR 分解），即

則X的r秩快速低秩逼近可改寫成

在實際操作中，常使用右投影矩陣Y1和左投影矩陣Y2來更新A1和A2，以有效提高低秩逼近的精度。式（7）中St取決于X-Lt的硬閾值，即

式中：?k(X-Lt)為矩陣元硬閾值算子，該算子保留|X- |Lt中最大的k個元素，而將其他元素置0。

綜上，雙邊隨機投影散斑去噪算法的求解步驟如下：

（1）輸入觀測到的OCT 數(shù)據(jù)X，設置參數(shù)q，tmax；

（2）初始化稀疏矩陣S、低秩矩陣L及高斯隨機矩陣A1，令L=X，S=0；

（4）計算右隨機投影矩陣Y1=L^A1，更新A2=Y1，計算左隨機投影矩陣A2（式（9））；

（5）更新A1，使得A1=Y2，若循環(huán)次數(shù)未達到q+ 1 次，則循環(huán)更新Y1和Y2，否則進行下一步；

（6）對Y1和Y2進行QR 分解（式（11）），求解Lt（式（12））；

（7）求解St（式（13））；

2 實驗結果與分析

實驗對象為10 名正常男性青年和5 名正常女性青年，對其以黃斑為中心的視網膜OCT 圖像數(shù)據(jù)進行測試。視網膜OCT 圖像均使用Topcon DRI OCT?1 掃描獲得；圖像大小為512 體素×992 體素×256 體素，覆蓋了以黃斑為中心的6 mm×6 mm×2.5 mm 區(qū)域。

2.1 圖像質量評價標準

為有效評價本文所提出方法的有效性，引入了常用的3 個圖像質量評價指標，分別為信噪比（Sig?nal to noise ratio，SNR）、對比度信噪比（Contrast to noise ratio，CNR）以及等效視數(shù)（Equivalent number of looks，ENL），其計算公式分別為

式中：MAX 為B 掃描圖像的最大可能像素值；σb和ub分別為背景區(qū)域方差和均值；σi和ui分別為若干感興趣區(qū)域的方差和均值；N為感興趣區(qū)域個數(shù)。SNR 值越高則代表OCT 圖像背景區(qū)域的圖像質量越好；CNR 反映了兩種組織信號強度的相對差別，主要用來衡量圖像前景和背景的對比度，數(shù)值越大則圖像的對比度越好；ENL 主要用來評價圖像光滑度和反映信號強度，與CNR 相似，ENL 通常選取若干感興趣區(qū)域進行評價。本文實驗在視網膜圖像神經纖維層、內視網膜以及視網膜色素上皮層上隨機選取了3 個區(qū)域作為感興趣區(qū)域。

2.2 仿真結果與分析

用本文所提的雙邊隨機投影去噪方法對臨床數(shù)據(jù)集進行測試，實驗結果如圖4 所示。在原始圖像和去噪圖像中選取對應位置的A 掃描圖像（圖4 中黃線位置）分別繪制一維深度信息信號圖（圖5）。隨著深度逐漸增大，視網膜從上而下對應的結構組織分別為神經纖維層（Retinal nerve fiber layer，RN?FL）、神經節(jié)細胞層（Ganglion cell layer，GCL）、內 從狀層（Inner plexiform layer，IPL）、內核層（Inner nuclear layer，INL）、外從狀層（Outer plexiform layer，OPL）、外核層（Outer nuclear layer，ONL）、內節(jié)/外節(jié)（Inner segments layer/ outer segments layer，IS/OS）、視網膜色素上皮層（Retinal pig?ment epithelium，RPE）以及脈絡膜（Cho?roid）等。對比圖5 中去噪處理前后的深度信息信號圖可見，去噪處理后信號曲線的光滑度以及視網膜生物結構信息對比度都有明顯提升。

圖4 雙邊隨機投影去噪實驗結果Fig.4 Experimental result of bilateral random projection denoising

圖5 一維深度信息信號圖Fig.5 One?dimensional depth information signal map

為更好地評價本文方法客觀性能，將原始圖像與本文提出的雙邊隨機投影去噪方法以及RPCA 算法［16］進行對比分析。 圖6為不同去噪方法處理所得的去噪效果對比圖。原始圖像以及各算法處理得到的去噪圖像的量化分析結果如表1 所示。

表1 各算法去噪結果Table 1 Denoising results of different algorithms

由圖6 可知，采用雙邊隨機投影去噪方法處理的圖像整體更加光滑、對比度更高、圖像質量更好。與RPCA 算法相比，本文方法在信噪比、對比度信噪比以及等效視數(shù)指標上分別提高了1.22 dB、0.84 dB 和59.5，去噪效果更好。這些指標的提高也意味著臨床醫(yī)生看到的OCT圖像中信號強度更強，視網膜各組織結構間對比度更好，各組織結構邊界更加平滑，圖像的可讀性更好。這將有利于提升醫(yī)生的閱片速度和閱片質量，同時也能有效提升計算機輔助診斷分析，如視網膜自動分層，病變區(qū)域自動分割等的精度。除此之外，本文對各方法的計算復雜度也進行了評估。 在Intel（R）Core（TM）2 Duo CPU 處理器和8 GB RAM 的PC 機，MATLAB 2017a 軟件上進行了測試，用RPCA 算法和本文算法對實驗數(shù)據(jù)集中單幅B掃描圖像的處理時間分別為4 s 和0.3 s。本文算法大大降低了計算復雜度。

圖6 各算法去噪結果Fig.6 Denoising results of different algorithms

為進一步說明本文方法的臨床價值，對10 名臨床早期青光眼患者的OCT 圖像進行了分割，并評估了散斑去噪對分割精度的影響。青光眼診斷的指標之一是神經節(jié)細胞層加上內叢狀層（GCL+IPL）以及神經節(jié)細胞復合層（RNFL+GCL+IPL，GCC）的厚度。人工對三維視網膜OCT 圖像進行分層測量費時費力，因此，眼科醫(yī)生需要一種可靠高效的方法來定量分析視網膜結構參數(shù)。為了評價散斑去噪處理對提升青光眼診斷精度和效率的作用，將去噪前后圖像分割的結果與分層金標準進行比較。金標準由兩名專業(yè)眼科醫(yī)生獨立手工標注，標注結果取平均值作為參考基準。由于三維OCT 圖像數(shù)據(jù)量很大，因此在標注時每個三維OCT 圖像中僅選取16 個B 掃描圖像進行標注。OCT 原始圖像，去噪后圖像均采用愛荷華大學提供的可視化軟件OCT?Explorer 進行自動層分割［21］。將分割結果與參考基準進行比較，計算層表面邊界定位誤差的均值，結果如表2 所示。結果表明，經本文算法去除散斑噪聲后，分割精度得到了提升。本文方法去噪后的圖像邊界定位誤差均值最小、分割精度最高，有助于醫(yī)生進行精準化的計算機輔助疾病診斷及分析。

表2 早期青光眼分割無符號表面誤差比較Table 2 Comparison of unsigned surface error in early glaucoma layer segmentation μm

3 結束語

針對含有噪聲的視網膜OCT 圖像，首先進行預處理，減少眼動偽差，從而保證相鄰幀生物組織結構之間的高度相似性；隨后，利用相鄰幀生物組織結構之間高度相似性及圖像高分辨率特性對原始OCT 圖像信號進行分解，將其分解為無噪低秩矩陣、稀疏矩陣以及噪聲矩陣；最后，采用雙邊隨機投影算法求解無噪低秩矩陣，從而達到去除散斑噪聲的目的。該方法的主要優(yōu)勢在于不需要大量訓練數(shù)據(jù)，利用少量相鄰幀即可在去除散斑的同時保留原有生物結構信息，運算效率較高。實驗結果表明，基于雙邊隨機投影的去噪算法能有效地抑制視網膜OCT 影像中的散斑噪聲，獲得了較好的視覺效果，且計算復雜度較低。除了應用在OCT 圖像去噪外，該方法還適用于高光譜影像去噪［22］、腦電信號提取［23］、雷達信號分析［24］和目標檢測等領域［25］，有廣泛的應用前景。該方法的提出為醫(yī)生進行視網膜疾病診斷以及改善計算機自動輔助視網膜疾病分析的性能提供了有效幫助，且計算時間可以滿足臨床的實時性要求。然而，對于病變視網膜，由于視網膜結構發(fā)生嚴重形變，相鄰幀的結構相關性也相應減少，從而影響算法精度，因此，在下一步研究工作中將重點考慮如何改進本方法使之適用于嚴重病變視網膜的去噪。