陸莉霞,鄒俊忠,郭玉成,張 見,王 蓓
1.華東理工大學 控制科學與工程學院,上海200237
2.清影醫療科技(深圳)有限公司,廣東 深圳518083
膝關節主要包括韌帶、關節面、肌肉、肌腱等多個組織結構,是人體結構中最復雜、最易發生損傷的關節之一。作為臨床上一種常見的骨損傷類型,膝關節損傷多發生在半月板、膝關節韌帶等位置,患者多表現為不同程度腫脹、疼痛等,且隨著現代交通業、建筑業的不斷發展,該疾病發生率呈現出逐年增高的趨勢,不僅給患者生活帶來了極大的不便,還影響到患者的生活質量[1]。
隨著X 線、CT、MRI 等影像學檢查的普及,膝關節損傷診斷的準確率顯著上升,為臨床早期治療提供客觀依據。MRI檢查具有較高的軟組織分辨率,能夠采集到更為豐富的圖像信息。通過矢狀位、橫軸位等方位的MRI 影像能夠對膝關節軟組織及周圍結構予以全面的探查,獲得詳細的、準確的圖像,因而MRI檢查對前十字韌帶撕裂、半月板撕裂的臨床診斷表現出較高的準確性[2-3]。近年來,膝關節磁共振成像(MRI)已成為診斷膝關節損傷的首選方法,一個用于分析膝關節MRI影像的自動化系統,可以篩選出高危患者,并幫助臨床醫生做出更準確的診斷。
傳統醫學影像分析主要采用邊緣檢測、紋理特征、形態學濾波以及構建形狀模型和模板匹配等方法。Haralick等人[4]使用灰度共生矩陣來提取乳腺組織的紋理特征;Swain 等人[5]使用顏色直方圖來提取眼底圖像的顏色特征。
MRI影像的多維度、多平面的特性限制了傳統醫學影像分析方法對膝關節MRI的應用[6],而深度學習方法能夠自動學習多層特征,非常適合于輔助診斷醫學影像[7]。近年來,深度學習方法已經超越了傳統的醫學影像分析方法,并在膝關節MRI 影像領域取得了重大進展。Prasoon等人[8]提出了一種新的體素分割系統,該系統集成了三個二維CNN 網絡,分別用于MRI 三維影像的xy、yz和zx三個平面的膝關節脛骨軟骨的分割。Liu 等人[9]開發了一套基于深度學習的全自動膝關節MRI軟骨損傷檢測系統,該系統有兩個CNN網絡組成,第一個CNN 網絡進行軟骨和骨骼的快速分割;第二個CNN網絡評估了分葉軟骨組織的結構異常。吳恩達團隊[10]首次提出了用于預測膝關節損傷的深度學習模型——MRNet,該模型通過MRI 影像對膝關節常見的損傷進行預測,提高MRI診斷的質量。
本文通過將膝關節矢狀面、冠狀面和軸向面的MRI影像的多種特征融合,對前交叉韌帶撕裂、半月板撕裂和一般異常三種常見的膝關節損傷進行分類預測。通過各個對比實驗,驗證本文方法的有效性和可行性。
對于膝關節MRI影像的特征提取方法,主要分為傳統特征提取方法和深度學習特征提取方法。一般來說,傳統方法提取的特征底層語義特征,例如邊緣特征和紋理特征。而深度學習如卷積神經網絡(Convolution Neural Network,CNN)等,提取的圖像特征包含更多的細節信息,能表達高層語義信息。底層特征語義信息與高層語義信息對于膝關節損傷的診斷均有重要意義,因此本文提出用于膝關節診斷的一種多模態特征融合網絡,如圖1所示。

圖1 多模態特征融合網絡流程圖
首先,提取邊緣特征HOG 特征和紋理特征LBP 特征,經contact 融合后利用PCA 選取特征貢獻度超過95%的特征作為傳統特征;其次,在VGG16網絡模型的基礎上加入金字塔融合的思想,將多個feature map 的信息融合作為深度特征;最后,通過一個多層神經網絡將傳統特征和深度特征經進行深度融合,并得到預測概率。
膝關節MRI 影像的一個樣本的數據大小為n×256×256,其中n為一組MRI影像序列所包含的靜態圖片數目。
由于膝關節主要位于MRI的中央位置,因此先將所有靜態圖片進行裁剪操作,并截取圖像中心224×224區域。接著進行歸一化操作,像素值被歸一化到[0,255]之間,歸一化是為了消除其他變換函數對圖像變換的影響。在訓練階段,數據進行數據增強,包括水平翻轉和旋轉角度為25°,水平和垂直方向移動0.1的仿射變換。
1.2.1 HOG特征
HOG特征提取算法是Dalal等[11]提出的一種特征提取算法,核心思想是通過計算圖像局部區域梯度,并將每個局部區域中各像素點梯度的方向直方圖級聯。圖2為HOG特征提取算法的基本流程圖。

圖2 HOG特征提取算法的基本程圖
HOG特征提取的具體步驟如下:
(1)為消除光照和噪聲的影像,首先對圖像進行圖像灰度化和伽馬標準化處理。
(2)用一維中心對稱算子k=[-1 01]及其轉置計算圖像橫坐標和縱坐標方向的梯度,表達式為:

其中,H(x,y)表示像素在(x,y)處的灰度值。得到像素點的梯度幅值和梯度方向為:

(3)將圖像分割為若干個小塊(block),每個小塊由4 個單元(cell)組成,每個單元由8×8 像素組成,塊的滑動步長為1個單元。將θ(x,y)在[0,π]區間上分為9個區間,即每20°為一個方向(bin)。單元(cell)中的每一個像素點都為直方圖通道進行加權投票,權重為g(x,y),這樣就得到每個單元內的9個方向的梯度直方圖。
最后按順序級聯所有block 的直方圖向量,得到圖像的HOG特征μHOG,維度為7 056維。
1.2.2 LBP特征
原始的LBP模式算子為3×3的局部鄰域,以鄰域中心點像素的灰度值I(c)為基準,對周圍其他8 個像素點做二值化處理,比I(c)大的鄰域像素值置為1,否則為0。將這8個二進制數按順序排列得到一個二進制序列即為中心點像素的LBP值。

為對原始的LBP模式進行降維,并在數據量減少的情況下能最好地表示圖像的信息,本文采用均勻局部二值模式(Uniform Local Binary Pattern,ULBP)[12]。當某個LBP 所對應的循環二進制數從0 到1 或從1 到0 最多有兩次跳變時,該LBP所對應的二進制就稱為一個等價模式,其余情況為非均勻模式。通過這樣的改進,模式數量由原來的2p種減少為p(p-1)+2,二進制模式的種類大大減少,而不會丟失任何信息。
本文用局部二進制編碼直方圖表示LBP 圖像的特征,將LBP 圖像分為8×8 塊區域,獲取每個區域的LBP編碼直方圖,繼而得到整幅圖像的LBP 編碼直方圖,一張LBP圖像共有3 776維特征,記為μLBP。
1.2.3 HOG特征與LBP特征融合
經上述方法提取到的HOG特征和LBP特征均存在特征維度過高,維度之間存在耦合的問題,本文通過主成分分析(PCA)[13]來達到降維的效果。PCA 不僅減少了特征集的維數,同時避免了重要信息的丟失,保持特征集中對方差貢獻最大的特征。
設向量X=(x1,x2,…,xn)T,PCA的具體步驟如下:
將HOG 特征和LBP 特征歸一化后,經contact 融合后的特征為[μHOG,μLBP]。經PCA 降維后,將特征向量按對應特征值大小從上到下按行排列,前1 325 項的累計貢獻率超過95%,為95.43%。因此將原特征集降維到1 325維得到傳統特征集μ=(μ1,μ2,…,μ1325)。
近年來,卷積神經網絡在醫療圖像診斷領域的應用得到了廣泛而迅速地發展,并已取得不錯的成果。卷積神經網絡通過對圖像的卷積操作,可自動完成特征的提取,且這些特征擁有更高級的語義信息,魯棒性更強[14]。
VGG16網絡[15]是2014年由牛津大學的視覺幾何組(Visual Geometry Group)提出,該網絡共有五個卷積塊,每個卷積塊后都是最大池化層,最后三層為全連接層,VGG16的結構配置見表1。

表1 VGG16結構配置
VGG16 網絡利用最后一層特征來進行分類,這種方法的優點是速度快、需要內存少。缺點是僅僅關注深層網絡中最后一層的特征,卻忽略了其他層的特征,但多種信息的融合可以在一定程度上提升診斷的準確率。由于卷積神經網絡中淺層可以學習到比較通用的低級特征,如基本的灰度和邊緣信息等,深層可以學到特征的更多細節。因此本文在VGG16 網絡的基礎上加入Top-down 結構,將提取的特征層逐步搭建為特征金字塔,并進行像素疊加,網絡結構如圖3所示。
通常深度學習被認為在大量有標注的圖像上進行訓練才可以達到良好的效果,在醫學圖像中很難有大量帶有標注的病理圖像用于分割訓練,因此本文采用遷移學習的方法,直接使用預訓練好的VGG16 的第二至第五個卷積塊卷積塊(Block)。將每個卷積塊的最后一個卷積層得到的特征映射提取出來,分別為F(n),其中n=2,3,4,5。F(n)經上采樣(upsampling)與F(n-1)經1×1 卷積處理后的結果進行像素點融合,再采用3×3 的卷積核對融合結果進行修正,消除了上采樣的混疊效應,得到新的特征映射F′(n-1)。金字塔融合的公式為:

獲得最后一層融合層的F′(2),依次通過BN 層、ReLU層、自適應最大池化層和全連接層。BN 層不僅能夠加快模型的收斂速度,還能提升分類效果[16]經過BN(Batch Normalization)層。設x(k)為F′(2)第k維特征,BN層為x(k)引入了兩個可學習的參數γ(k)和β(k),并利用無差估計輸出第k維特征:

ReLU 層加入非線性因素的,增加模型的表達能力不夠。ReLU激活函數的表達式為:

自適應最大池化層與標準的Max Pooling 區別在于,自適應最大池化層Adaptive Max Pooling 會根據input_size(In)來控制輸出output_size(Out)。stride和κernelsize分別為:

全連接層可以看做h×w的全尺寸卷積,h和w為前一層輸出的大小,最終得到1 024 維通過卷積神經網絡提取的深度特征g=(g1,g2,…,g1024)。
由于提取到的傳統特征和深度特征維度均較高,若直接級聯兩個特征向量作為最后的特征會大大增加時間成本和計算成本,因而在多模態特征融合時要減少融合后特征的維度。同時考慮到不同模態的特征在融合結果中所起到的作用,以及保留多模態的相關性,本節構建了一個多模態特征融合模塊,該模塊包含一個神經元個數小于特征維數的隱藏層和一個Sigmoid 層,通過最大化特征層的能量比重來訓練整個網絡,如圖4 所示。隱藏層對融合后特征進行非線性降維,Sigmoid 層將融合后特征映射到(0,1)區間,即預測概率。令特征向量x=(μ,g),前向傳播公式為:


圖3 加入金字塔融合的VGG16網絡結構圖

圖4 多模態特征融合網絡圖

其中,T為樣本數量。當函數ψ(θ)取得最大值時,特征層的能量比重大,隱層的能量小。當在網絡內傳輸數據時,數據流的方向也是能量衰減的方向,多次迭代之后,網絡能量呈衰減趨勢,網絡趨于有序或者概率分布趨于集中。
本文將偏差βi、αj初始值設為服從平均值為0,標準差為0.1的正態分布的隨機數,權重γij初始值設為服從平均值為0,標準差為0.01的正態分布很小的隨機數,隱藏層神經元設為1 000個。
實驗時,所使用的CPU 型號為3.4 GHz、16 核AMD1950x,內存容量為64 GB,GPU 為12 GHz 顯存NVIDIA TITAN-V。軟件方面,使用Python作為設計語言,配合opencv等視覺庫進行代碼編寫。深度學習模型框架使用Pytorch以完成對模型的構建,訓練與預測。
本文網絡模型使用Adam優化器,增加模型收斂穩定性和收斂速度,同時也能得到優良的分類結果。學習率經過調試設置為1E-5,權重衰減為0.01,模型運行50輪。
實驗過程中使用了公開的膝關節MRI 競賽數據集[10]對本研究中的網絡模型進行訓練和驗證。本研究的核心任務為檢陽,因而本文采用了以下評價指標:準確率Accuracy、Recall 和AUC(Area Under Curve)值定量的評估不同模型的性能。
準確率為模型預測正確的樣本占所有參與預測的樣本的比例。定義為:

其中,TP為真陽性樣本數量,TN為真陰性樣本數量,FP為假陽性樣本數量,FN為假陰性樣本數量。
AUC 值被定義為ROC 曲線下與坐標軸圍成的面積。由于ROC曲線一般都處于y=x這條直線的上方,所以AUC 的取值范圍在0.5 和1 之間。AUC 越接近1.0,檢測方法真實性越高。
本次研究中的膝關節MRI 影像數據集來自于斯坦福大學醫學中心整理的從2001 年1 月1 日至2012 年12 月31 日期間內的1 370 個膝關節核磁共振檢查數據集——MRNet數據集[18]。該數據集包含1 104例(80.6%)異常(Abnormal)檢查,其中前交叉韌帶撕裂(ACL tear)319例(23.3%),半月板撕裂(Meniscus tear)508例(37.1%)。194例(38.2%)患者同時發生前交叉韌帶撕裂和半月板撕裂。每個樣本包括矢狀面加權序列(sagittal plane)、冠狀面加權序列(coronal plane)和軸向面加權序列(axial plane)。這些序列的圖像數量從17 幅到61 幅不等(平均31.48幅,SD7.97幅)。
數據集被分為訓練集(1 130個樣本,1 088名患者)和驗證集(120個樣本,111名患者),訓練集與測試集的數據分布見表2。

表2 訓練集與測試集的數據分布
若將矢狀面加權序列、冠狀面加權序列和軸向面加權序列同時放入上述多模態特征融合網絡中,不僅對硬件要求較高,運行時間較長,而且不同序列的特征表達差異較大,預測準確度不高。
因而本研究針對每個任務(異常、前交叉韌帶撕裂、半月板撕裂)和序列類型(矢狀面、冠狀面、軸向面)訓練了不同的9個模型。
在驗證階段,采用回歸模型得到異常、前交叉韌帶撕裂、半月板撕裂的預測概率,回歸模型結構圖如圖5所示。

圖5 回歸模型結構圖
將同一個樣本的矢狀面、冠狀面和軸向面加權序列同時放入上述多模態特征融合網絡中,分別得到在該序列下,樣本為異常樣本的概率P(A)、P(B)和P(C)。然后,將P(A)、P(B)和P(C)分別乘以一個回歸系數,相加后通過Sigmoid 函數得到在三種序列下,此樣本為異常樣本的預測概率,計算公式為:

為了尋找最佳參數wA、wB和wC,本文采用梯度上升的方法,將梯度算子移動的步長記為α,則梯度算法的迭代公式為:

同樣地,可以獲得此樣本為前交叉韌帶撕裂樣本和半月板撕裂樣本的預測概率。
特征提取不僅僅是分類和預測的前提,也是分類準確率的保證,而高效準確的特征可以提高后續分類和預測的精確度。本文數據經相同的預處理后,每種膝關節MRI序列提取的傳統特征(Model a)、CNN特征(Model b)和多模態融合特征在每個任務下的實驗結果見表3。
由表3 中Model a 的結果可知,僅僅依靠MRI 影像的紋理和邊緣特征,在膝關節損傷預測任務中存在一定的局限性,尤其在半月板撕裂的預測上,最高的準確率僅有69.17%,無法很好地輔助臨床醫生進行預測。而與傳統方法相比,深度學習方法在膝關節MRI影像診斷領域具有明顯優勢。由表3中Model a和Model b的結果可知,深度學習方法的準確率、召回率和AUC值均有大幅提升,其中在異常預測任務中召回率最高已達到97.89%,同時在傳統方法有所欠缺的半月板撕裂預測任務中,深度學習方法的準確率均高于76%,充分顯示了深度學習方法能夠更好地診斷出陽性的樣本,同時AUC值也更高,因而在偏態的樣本中更穩健。

表3 各模型效果對比
本文模型在卷積神經網絡中融入低級語義信息,豐富了模型的特征表達。由表3中Ours的結果可知,在融合了傳統特征與深度學習特征后,各性能指標均得到一定程度的提升,尤其是準確率和AUC值。在異常、前交叉韌帶撕裂、半月板撕裂這三項任務中,最高準確率為92.50%、92.50%和81.83%,最高AUC 值為0.961 3、0.974 2 和0.857 4。同時,本文模型對一般異常的陽性樣本的檢出率均高于94%,對前交叉韌帶撕裂和半月板撕裂陽性樣本的檢出率也較高。這對于優先發現高危患者并協助臨床醫生進行診斷有著重要意義。
由表3 中Ours 的結果可知,在一般異常預測任務中,軸向面序列表現最佳,準確率、召回率和AUC 值分別為92.50%、97.95%和0.961 3;在前交叉韌帶撕裂預測任務中,矢狀面序列更有優勢,準確率、召回率和AUC值分別為92.50%、94.44%和0.972 8;在半月板撕裂預測任務中,軸向面序列表現較好,準確率、召回率和AUC值分別為81.83%、82.92%和0.846 5。此結果與半月板位于兩側脛骨與股骨外側髁間[19]而前交叉韌帶位于膝關節滑膜外側[20]的膝關節結構相符合,因此,本文實驗結果具有一定的科學性。
鑒于每個任務均只有陽性和陰性兩種可能,為了更好地展示回歸模型的效果,將回歸模型預測概率大于或等于0.5 的認為是陽性樣本,小于0.5 的認為是陰性樣本,則在一般異常、前交叉韌帶撕裂和半月板撕裂下的實驗結果如表4。

表4 回歸模型各性能指標
由表4 可知,本文模型在一般異常、前交叉韌帶撕裂和半月板撕裂預測任務中均有較好的性能指標。尤其是一般異常和前交叉韌帶撕裂預測任務,準確率和召回率均超過90%,AUC 值均高于0.94,說明本次研究所應用的方法能夠在這兩項膝關節損傷預測中為臨床醫生提供可靠的基礎診斷。由于半月板主要由纖維軟骨構成,而MRI影像對軟骨組成成分具有局限性[21]。在這種情況下,本文模型在半月板撕裂預測任務中的AUC值仍達到了0.847 9。綜合來看,本文方法在膝關節損傷診斷方面具有一定的價值。

表5 不同膝關節損傷檢測模型性能比較
表5 對比了本文算法與MRNet 競賽[18]前三名結果和現有文獻中最好結果的平均AUC值。
可以發現,在對膝關節損傷進行預測時,表5中6種深度學習網絡AUC 值均高于90%,都能夠取得較好的預測結果,驗證了深度學習強大的特征學習能力。Bien等[10]在論文中提出的MRnet 以及mrnet-baseline 都是以alexnet 模型為基礎的,AUC 值雖已達到0.917,但受alexnet模型本身深度和廣度的限制[22],卷積層能學習到的特征的數量和種類有限。本文模型以更深層的vgg16模型為基礎,并在該卷積神經網絡的基礎上加入了多模態的融合特征,包括多層卷積層的融合特征以及MRI影像的紋理和邊緣特征,讓模型有了更強大的表達能力,以此來獲得更好的預測效果。因此本文模型的平均AUC值略高與其他深度學習網絡,為0.920。
本文提出了用于膝關節診斷的一種多模態特征融合網絡。首先,通過對膝關節MRI影像的預處理,從傳統與深度學習兩方面提取膝關節損傷的多模態特征;然后,通過多層神經網絡對特征進行快速的相關性融合;最后,用logistic 回歸整合各個模型的概率。通過不同的實驗,對上述所提到的三種模型進行了比較,同時對回歸模型進行了分析,結果顯示本文方法在膝關節損傷預測方面具有較大的優勢和較高的實用價值。在未來的研究工作中,將嘗試更多模態的特征融合,例如MRI影像在時間維度上的特征。同時將此應用在膝關節MRI 的算法應用在其他醫學圖像類型和不同人體器官的診斷上。