999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于數字全息的血紅細胞顯微成像技術*

2020-08-29 07:28:24張益溢吳佳琛郝然金尚忠曹良才
物理學報 2020年16期

張益溢 吳佳琛 郝然 金尚忠? 曹良才?

1) (中國計量大學光學與電子科技學院, 杭州 310018)

2) (浙江省現代計量測試技術及儀器重點實驗室, 杭州 310018)

3) (精密測試技術及儀器國家重點實驗室, 清華大學精密儀器系, 北京 100084)

1 引 言

數字全息顯微技術(digital holographic microscopy, DHM)結合數字全息與顯微成像技術, 充分利用兩者優勢, 不僅能夠對待測樣本進行放大, 而且能實時記錄待測樣本的三維信息.1999 年 E. Cuche 課題組在待測樣本與CCD(charge coupled device)像感器之間放置顯微物鏡[1], 用于放大待測樣本的物光波, 僅基于單幅全息圖便可獲取有關樣品三維結構的精確定量信息, 降低了對CCD 分辨率的要求, 實現了高分辨率顯微成像. 但由于顯微物鏡的小孔徑限制, 全息圖部分高頻成分在記錄過程中丟失, 重建像質量受到影響, 因此數字全息顯微技術仍存在需要改進的地方. 2003 年Ferraro 等[2,3]研究了自動跟蹤聚焦的數字全息顯微鏡, 并在隨后的研究中提出了運用多波長數字全息技術來校正色差. 2006 年Charrière 等[4]將數字全息顯微技術與電子計算機斷層掃描技術(computed tomography, CT)相結合, 通過旋轉花粉細胞獲取其各方向的相位圖分布, 根據濾波反投影算法實現了生物樣品的三維成像. 2007 年Kemper 等[5]將數字全息顯微技術與生物醫學相結合, 重建了人體肝臟腫瘤細胞與血紅細胞(red blood cell, RBC)的三維形貌, 進一步拓展了數字全息顯微的實踐應用. 2008 年Kim 等[6]通過數字全息顯微術結合生物醫學, 揭示了血管的形貌以及視網膜層的光學厚度分布. 國內圍繞數字全息成像分辨率提升與生物細胞定量相位測量開展了大量的工作. 在分辨率提升方面, 2010 年翟宏琛等[7]利用一種相位掩模板進行多角度照明,實現了綜合孔徑數字全息超分辨率成像;2013年馬俊等[8]基于結構光照明的數字全息顯微技術,通過分離和合成傅里葉域中不同頻率區域的帶寬,可以獲得更高空間分辨率的重建圖像;同年王華英等[9,10]綜合分析了數字全息顯微成像系統的分辨率影響因素,對像面數字全息系統展開了討論與研究;2015年姚保利等[11]利用散斑照明獲得重建物波中的合成數值孔徑,以提高分辨率;同年王大勇等[12,13]利用動態空間光柵提高了數字全息成像分辨率,并通過微球體與像面數字全息相結合,實現了分辨率增強的相襯成像.在上述研究中,由于分辨率得到了大幅提升,使得數字全息顯微更有效地分辨生物微觀結構如人體細胞與亞細胞結構等.在數字全息定量相位測量方面,2007 年董可平等[14]成功實現了生物細胞的定量相位測量;2008年趙建林等[15]通過拍攝兩幅全息圖并利用相位相減技術實現了活體細胞組織等相位型生物樣本的定量測量和有效觀察;2012年馬利紅等[16]建立了一套預放大式數字全息顯微成像系統,實現了高分辨率的細胞定量相位測量;2018年肖文等[17]利用離軸馬赫澤德實驗系統對細胞進行實時成像,證明了細胞質中肌動蛋白對細胞骨架具有更大的支撐作用.2019年左超等[18,19]通過數字全息顯微結合自適應松弛超像素分辨技術,在利用相位恢復與數字重聚焦等算法或方法的基礎上,開發了無透鏡數字全息顯微系統,集成了定量相位、剖面分析、三維顯示和細胞計數等多項功能,具有無標記、高通量、低成本和微型化的特點,為面向臨床檢驗應用的顯微設備提供了有效的技術支撐.基于數字全息的定量相位成像技術滿足了活體細胞組織等相位型生物樣本定量測量和動態觀察的需求,對于生物醫學領域的發展具有重要意義.

相較于當前細胞成像領域廣泛應用的激光掃描共聚焦顯微(confocal laser scanning microscope,CLSM)技術[20]、近場光學顯微(near-fieldoptical microscopy,NSOM)技術[21]和光學相干層析成像(opticalcoherencetomography,OCT)技術[22]等,

數字全息顯微技術具備以下優點:(1)屬于面成像,在成像速度上極具優勢,僅根據單次曝光即可實時記錄待測樣本完整三維信息[5,23],利于生物活體組織的動態變化的研究,傳統方法主要基于點掃描成像;(2)屬于寬光束照射,對光功率密度要求較低的同時無需熒光標記,實現了無擾分析,而激光掃描共聚焦需在樣品上實現激光束聚焦,對焦點處的光功率密度要求較高,或需要對樣品進行熒光標記,無法實現無擾分析;(3)可以得到細胞折射率分布信息,對比當前CT技術具有更高的空間分辨率,且對比度高,可實現自動聚焦,近場光學顯微技術僅能實現表面成像,而光學相干層析成像技術雖能得到生物體組織內部的折射率分布信息,其空間分辨率約為5μm,成像深度在毫米量級,需要進一步改進才能觀察人體細胞和亞細胞結構(1μm左右).通過數字全息顯微技術能夠實現寬視場亞微米級空間分辨率(0.5μm 左右),不僅可以清晰分辨人體細胞與亞細胞結構,并且可實現定量、三維和快速追蹤的成像,在細胞成像領域上具有獨特的應用價值.

本文從生物醫學細胞成像應用需求出發,介紹了數字全息顯微成像的原理,綜述了同軸、離軸以及光鑷輔助離軸數字全息顯微系統如何對RBC進行形貌信息提取,重點分析了瑞利索末菲反向傳播算法、清晰度量化算法、分水嶺分割算法、數字重聚焦方法與熱漲落方法在研究RBC的微形變、三維體積測量與空間分布等方面的應用,概述了當前數字全息顯微成像的發展與應用前景.

2 數字全息顯微成像原理

2.1 數字全息顯微記錄原理

數字全息顯微光路如圖1所示,其中Object plane 代表物平面x0-y0,MOplane 代表顯微物鏡平面xm-ym,Hologramplane代表全息記錄平面x-y,Imageplane 代表透鏡成像平面xi-yi,樣品放置于距離顯微物鏡平面z0距離處[24],顯微物鏡平面與全息記錄面距離為zm,全息記錄面與透鏡成像平面的距離為z,R 作為參考光,與經顯微物鏡放大后的物光O于全息記錄面上發生干涉,形成明暗相間的干涉條紋,這些條紋記錄了物體光波所有振幅與相位信息,是一個編碼和調制的過程.

設物光在物平面的光場分布為O0(x0,y0),根據菲涅爾衍射近似條件,物光波經過距離z0衍射至顯微物鏡前表面時, 其光場復振幅分布Of(xm,ym)為

圖1 數字全息原理示意圖Fig. 1. Optical layout of digital holography.

其中,λ表示物光波長,x0與y0表示物平面上的空間坐標,xm與ym表示顯微物鏡平面上的空間坐標,z0表示物平面與顯微物鏡平面的距離. 由于顯微物鏡在全息記錄過程中僅作為放大作用的透鏡, 依據透鏡的相位調制特性, 物光波經過顯微物鏡傳播至其后表面時, 其光場復振幅分布Ob(xm,ym) 為

其中,f為顯微物鏡的焦距, 物光波由顯微物鏡后表面經過距離zm傳播至全息記錄面時, 其光場復振幅分布O(x,y) 為

由于顯微物鏡的放大作用, 全息記錄面記錄的物光場相對原物光場存在坐標縮放的關系. 同理,參考光傳播至全息記錄面的光學路徑與物光類似,參考光到全息記錄面的距離與物光到全息記錄面的距離相等, 設參考光在全息記錄面的光場分布為R(x,y) , 物光和參考光在全息記錄面上發生干涉形成全息圖, 全息圖的光場分布為

則全息記錄面上全息圖的強度分布為其光場分布的平方, 其強度分布為

式中, 第一項表示物光的強度分布, 第二項表示參考光的強度分布, 第一項與第二項合稱為零級項或直流項; 第三項表示物光場的復振幅分布, 包含了物光波波前的振幅信息與相位信息, 進行數值重建時可獲得待測樣本的真實物像[25]; 第四項表示物光場的的共軛復振幅分布, 數值重建時會產生孿生像, 這是全息圖重建時獲得重建像的最大干擾項,需要在重建過程中通過算法進行分離.

2.2 數字全息顯微重建原理

數字全息重建通過計算機數值模擬重建光波經過數字全息圖后的衍射過程, 實現數值重建, 在數字全息顯微成像系統中重建的并非物光本身, 而是由顯微物鏡放大后所成的物光場, 根據近似條件不同可將模擬衍射成像方法分為菲涅爾重建法、卷積重建法、角譜重建法[26]. 三種重建法均基于菲涅爾-基爾霍夫衍射公式進行近似處理或變換后所得.菲涅爾重建法重建圖像的像素尺寸δx與重建距離z滿足關系:δx=λz/X0, 其中X0為全息圖尺寸.這意味著重建圖像的大小會隨著重建距離增加, 而采用卷積重建法和角譜重建法所得圖像像素尺寸始終與全息圖像素尺寸相同. Kreis 等[27]采用菲涅爾重建法與卷積重建法分別對CCD 記錄下的骰子全息圖進行了重建, 其中骰子位于CCD 前1.054 m處, 且骰子尺寸大于CCD 尺寸. 采用菲涅爾重建法可直接重建完整的骰子圖像, 而采用卷積重建法需要平移CCD 采集多幅全息圖或者補零操作才能完整重建骰子圖像, 因此菲涅爾重建法適用于尺寸大于像感器的待測樣本. 王大勇等[28]分別利用菲涅爾重建法、卷積重建法與角譜重建法對骰子全息圖進行重建, 比較了三種方法的重建速度, 其中菲涅爾重建法的執行時間為0.97 s, 速度最快, 而卷積重建法的執行時間為2.36 s, 速度最慢, 另外角譜重建法的執行時間為1.22 s. 菲涅爾重建法僅通過一次正傅里葉變換即獲得重建像; 卷積重建法由于其傳播函數在空域中表示, 需要經過兩次正傅里葉變換和一次逆傅里葉變換獲得重建像; 角譜重建法與卷積重建法原理類似, 其傳播函數在頻域中表示, 需要經過一次正傅里葉變換和一次逆傅里葉變換獲得重建像, 通過比較像質得到角譜重建法的準確性優于卷積重建法. 而根據USAF 鑒別率板的重建, 得到在卷積重建與角譜重建中均存在最優重建距離, 當重建距離不同于最優距離時, 重建像的分辨率會下降, 此時角譜重建法的結果優于卷積重建法[29,30]. 卷積重建法僅能在最佳重建距離附近較小的范圍內獲得高分辨率重建像, 適用性上存在不足, 而角譜重建法重建速度較快, 在最佳重建距離附近較大范圍都能獲得高分辨率的重建像[31].

下面以卷積重建法與角譜重建法為例, 介紹生物細胞全息圖的重建, 卷積重建過程由記錄全息圖的光場分布H(x,y) 以及卷積重建空域內的傳播函數h(xi ?x,yi ?y) 構成, 通過快速傅里葉變換(fast Fourier transform, FFT)進行運算:

而角譜重建過程由光場分布H(x,y) 與角譜重建頻域內的傳播函數G(fxi,fyi) 構成, 通過FFT 進行運算:

其中,F與F?1分別表示快速傅里葉變換與快速傅里葉逆變換, 如圖2 所示,x,y表示全息記錄面上的空間坐標,xi與yi表示重建成像平面上的空間坐標,fxi與fyi表示成像面上 (xi,yi) 對應的頻率坐標,K(fx,fy)=F{H(x,y)}表示在全息記錄面上的全息圖頻譜.

圖2 數字全息顯微重建原理圖Fig. 2. Optical layout of digital holographic reconstruction.

重建圖像的強度可根據|Γ(xi,yi)|獲得, 由角譜分析方法可對(8)式進行進一步簡化[32], 其表達式如下:

其中,λ表示光源波長,z表示全息記錄面與成像面的距離. 可通過(7)式對記錄面所得全息圖進行卷積重建或結合(8)式與(9)式對記錄面所得全息圖進行角譜重建, 獲得經顯微物鏡放大后待測樣本的數值重建像.

對于同軸全息系統, 由于重建像與孿生像重疊, 需要在數字全息記錄過程中引入相移裝置, 再根據相移法將所需重建像與孿生像進行分離, 提取出待測樣本的物像[33,34]; 對于離軸全息系統, 選擇適當的物光與參考光的入射夾角, 使物光與參考光以一定夾角入射至像感器記錄面上, 產生干涉全息圖, 通過濾波以濾除孿生像, 提取出待測樣本的物像. 在離軸全息系統中, 由于物光與參考光存在夾角, 會產生一次項附加相位(呈現斜面分布), 而無論是離軸還是同軸全息系統, 當RBC 等弱散射物體的相位信息加載于經顯微物鏡后形成的球面波上時, 均會引入二次項附加相位(呈現球面分布),這些附加相位所產生的調制會使待測樣本自身的相位分布難以恢復, 引起包裹相位的欠采樣, 最終影響到相位解包裹的準確性, 可在參考光路放置與物光光路完全相同的顯微物鏡進行補償[35], 或通過數據擬合進行自動相位校正[36], 或采用前后兩次拍攝全息圖進行相位相減來實現相位校正[15],此時通過同軸或離軸系統獲得的校正相位分布在[–π, π]之間, 還需要進行相位解包裹以實現待測樣本相位圖像的恢復[37,38].

為有效測量RBC 的微形變, 三維體積及其空間分布, 當前相關RBC 測量的數字全息顯微系統在重建時主要采用卷積重建法或角譜重建法.

3 血紅細胞數字全息顯微成像系統

RBC 是血液中的主要有形成分, 約占血液體積的50%, 生存環境的改變會引起RBC 病理性變化, 如糖尿病患者體內的高血糖環境會導致RBC 雙凹圓盤狀喪失、形態拉長、體積減小、變形性降低[39,40], 心血管疾病患者體內RBC 體積分布寬度(RDW)升高[41,42], 帕金森氏疾病患者體內存在氧化應激, 氧化應激會引起體內RBC 功能性障礙與凋亡[43,44], 因此如何提取RBC 的生物學參數與形貌信息將有助于研究人體的糖尿病、心血管疾病、帕金森氏疾病等. 常規的RBC 血液測量技術血管造影與超聲多普勒分析所提供的血管形狀信息與血流速度信息的空間分辨率較差[45,46]; 粒子圖像測速(particle image velocimetry, PIV)在血液動力學信息的測量得到廣泛應用, 但其僅能在較淺的景深中固有地提供二維平面信息, 在應用范圍上存在較大缺陷[47,48]; 光學相干層析成像在血流動力學領域中得到廣泛應用, 能實現生物組織微觀結構顯示,但RBC 的對比度較低[49]; 光聲成像(photoacoustic imaging, PAI)能同時實現灌注成像與血流成像,且RBC 的對比度較高, 但其難以完成單個RBC 的動態追蹤[50]. 本節綜述了同軸、離軸以及光鑷輔助離軸數字全息顯微技術, 這些技術利用瑞利索末菲反向傳播算法、清晰度量化算法、分水嶺分割算法、數字重聚焦方法等不僅可以獲取高分辨率的血管形狀與血流速度信息[51], 進一步地還能夠檢測單個RBC 的形狀變化[52], 甚至可以完成RBC 動態追蹤, 實時測量RBC 三維體積, 實現無損3D 成像[53,54], 對于人體的糖尿病、心血管疾病、帕金森氏疾病等病理研究具有重大意義.

3.1 同軸數字全息顯微系統測量RBC 形變

為實現以數字全息顯微成像技術為基礎, 高效率、低計算強度、低成本地分析細胞的形態變化并跟蹤軸向位置提取三維坐標, 瑞典于默奧大學Andersson 等[52]設計了一套同軸數字全息系統并提出了一種方法, 通過瑞利索末菲反向傳播算法重建復振幅實部信息, 精確確定RBC 軸向位置并提供幾何形狀信息. 其光學系統的分析與實現可根據圖3 所述結構進行描述, 光學系統圍繞Olympus IX71 顯微鏡引入相關光學元件, 采用低成本發光二極管(470 nm, M470 L3-C1)發射LED 光[55], 多模光纖安裝于壓電控制顯微鏡載物臺(PS)上方[56],載物臺上放置單個RBC 樣品, 其通過人體血液與500 μL 的 pH=7.4 磷酸鹽緩沖溶液混合并離心提取制備而成, 多模光纖可以對LED 光進行頻率濾波, 穩定輸出光源, 經過RBC 發生衍射的光作為物光, 而未經過待測樣本的光作為參考光, 兩束光均 經 過 浸 水物鏡MO(60×/1.40NA)放 大,在CCD 像感器上發生相干疊加, 記錄下RBC 全息圖, 其中顯微鏡的壓電控制載物臺可進行軸向移動以精確控制物面與記錄面的距離 z , 以10 μm 為記錄間隔分別記錄z = 30—100 μm 共8 幅RBC 全息圖像. 根據同軸全息重建算法分別對8 幅RBC 全息圖進行重建, 計算各RBC 全息圖像在不同重建距離 z′下的歸一化強度, 認定歸一化強度最大位置即幾何焦點位置, 如圖4(a)所示, 插圖表示仿真實驗使用6.3 μm 寬, 折射率 n=1.40 的RBC形狀的Cassini 模型, 綠線與紅線分別代表實際實驗與仿真實驗中RBC 在不同記錄距離z =30—100 μm下對應的最佳反向重建距離, 可以看出, 在不同的記錄距離z 下, 通過分別計算RBC全息重建圖的歸一化強度最大值, 可以精確地得到8 個記錄距離下對應的幾何焦點位置, 作為最佳反向重建距離, 仿真實驗與實際實驗所得到的反向重建距離曲線一致.

圖3 同軸數字全息顯微成像系統檢測血紅細胞形變[52]Fig. 3. Inline digital holographic microscopic system for detecting micro-deformation of RBC[52].

第一步, 在合適的記錄距離下, 將實際實驗RBC 與仿真實驗RBC 在不同重建距離z′ = 0—100 μm 復振幅實部信息 Re(U) 進行對應比較, 圖4(b)中表示綠線與紅線分別代表當記錄距離z = 60 μm時, 實際實驗RBC 與仿真實驗RBC 的復振幅實部信息 Re(U) 隨反向重建距離變化而相應的值, 在相變位置z′ = 60 μm 周圍, 實際實驗RBC 與仿真實驗RBC 的 Re(U) 曲線吻合良好; 第二步, 進一步研究RBC 微小形態變化對 Re(U) 的影響, 通過改變仿真實驗RBC 的模型參數, 使初始RBC 發生微小形變, 微形變后RBC 在不同記錄距離 z 下歸一化強度峰值分布與初始RBC 基本重疊, 但其Re(U) 變化較大, 如圖4(b)中紫線所示, z′ = 60 μm附近波峰值提高, 而波谷位置由z′ = 50 μm 處移至z′ = 35 μm 附近, 且波谷值出現了下降, 表明RBC 微小形變可根據檢測 Re(U) 的變化得出; 第三步, 為證明提出實部信息 Re(U) 重建方法的魯棒性, 作者在實際實驗中選取了比初始RBC 體積小約20%的新型RBC 作為對象, 在仿真實驗中更改RBC 模型參數, 使其表示為一個5.04 μm 寬的RBC 形狀模型, 其實際實驗RBC 與仿真實驗RBC 在不同重建距離z′ = 0—100 μm 復振幅實部信息 Re(U) 的曲線如圖4(c)所示, 實驗結果表明, 在選取不同體積的RBC 的情況下, 實際實驗與仿真實驗的 Re(U)曲線仍良好吻合, 因此 Re(U)提供了足夠信息以區分細胞不同的幾何形狀. 綜上所述, 所提出瑞利索末菲反向傳播算法可以實現細胞形狀的辨別, 應用于生物學中數量龐大且深度不同的細胞在特定時間段內形態變化的病理研究.

圖4 RBC 形變與重建的復振幅實部信息之間聯系的結果分析 (a)不同記錄距離z 下最佳反向重建距離z', 綠線代表實際實驗, 紅線代表模擬實驗, 插圖代表模擬實驗所用RBC; (b)綠線代表實際實驗RBC 的Re(U), 紅線與紫線分別代表模擬實驗中初始RBC 與變形后RBC 的Re(U); (c)所用RBC 比(b)所用小約20%, 綠線與紅線分別代表實際實驗與模擬實驗所得Re(U)結果;(b)與(c)所用黑線表示RBC 的記錄距離[52]Fig. 4. The relationship between the RBC deformation and the real part of reconstructed amplitude. (a) The reconstruction distance z' to the focus at different recording distances z of the simulated RBC (red) and the experiment using a real RBC (green). Inset shows the Cassini model of the RBC used in the simulations; (b) the experiment using a real RBC (green), the reconstructed Re(U) for the simulated unaltered RBC and a deformed RBC represented by the red and blue curves, respectively; (c) reconstructed data from a ~20% smaller RBC compared with the one used in (b). Red and green curves represent the Re(U) of simulation and experiment, respectively; Gray vertical line in (b) and (c) indicates position of the RBC[52].

3.2 同軸數字全息顯微測量RBC 空間分布

通過高空間分辨率地提取血流量信息, 進而實現人體循環系統血管疾病的跟蹤監測, 這是當前生物醫學領域的重要課題之一. 韓國浦項科技大學Lee 等[51]提出了一套同軸數字全息顯微系統, 通過數字全息顯微成像技術對流動的RBC 逐一進行追蹤, 引入了清晰度量化算法對重建的RBC 圖像進行評估, 定位RBC 的精確深度位置, 可獲取單個RBC 的運動狀態, 實現高空間分辨率的血流量信息測量. 其光學系統的分析與實現可根據圖5 所述結構進行描述, 注射泵(SP)注射血溶液經過肝素化處理, 防止凝結, 溶液中RBC 數量為3.2 ×103cell/mm3以防止RBC 在全息圖中出現重疊.He-Ne 激光器發出632.8 nm 波長紅光, 經過水溶液與氟化乙烯丙烯(FEP)微管后進入血溶液, 部分光經過血溶液中的RBC 發生衍射作為物光, 而未通過RBC 而直接透過的光作為參考光, 兩束光由浸水顯微物鏡MO(20×/0.50NA)放大后在CMOS 記錄面上發生相干疊加, 產生干涉全息圖,由于血溶液折射率nmedium=1.346 , FEP 微管折射率nFEP=1.338, 水的折射率nwater=1.333 , 三者折射率幾乎相同, 避免了折射像畸變問題. 實驗初始階段首先移動CMOS 像感器以找到同時滿足成像清晰與足夠視場兩個條件的合適位置, 將CMOS像感器固定, 在測量體積600 μm × 350 μm ×350 μm 范圍內可觀察到185 個RBC. 根據同軸全息重建算法對RBC 圖像進行重建, 重建距離設定為深度方向上以2 μm 為間隔建立共計200 個重建切片進行RBC 圖像的重建得到重建圖像.

圖5 同軸數字全息顯微成像系統追蹤紅細胞空間分布[51]Fig. 5. Inline digital holographic microscopic system for tracking spatial distribution of RBCS[51].

為得到流動RBC 實際深度位置z, 引入圖像清晰度量化算法對重建圖像進行判斷[57]:

式中表示重建圖像強度梯度的變化狀況, 其中參數指標 LAP(z) 可以在忽略恒定強度梯度部分的情況下, 用于測量圖像的高頻邊緣成分進而反映RBC 的聚焦值, 在本實驗中對比其余圖像清晰度量化算法其精度更高[58,59], 參數指標 LAP(z) 的大小隨著重建距離z′的變化而變化, 可用于描述同一個RBC 在不同距離的重建切片下其聚焦值的大小變化, 由于實驗在深度方向上以2 μm 為間隔建立共計200 個重建切片, 則同一個RBC 在200 個重建距離下存在200 個聚焦值, 以重建距離為橫坐標, 用光滑曲線表示200 個聚焦值, 其聚焦值最大位置表示焦點位置, 定義此位置為RBC 的最佳重建距離, 根據計算顯微物鏡的軸向放大率可進一步得到RBC 的深度位置. 圖6(a)表示CMOS 像感器所采集的RBC 全息圖, 圖中共存在185 個RBC, 且RBC 的深度互有異同. 首先對RBC 進行分割與重建, 根據清晰度量化算法對重建圖像進行最佳重建距離判定, 進一步獲得RBC 的深度位置. 圖6(b)—6(d)分別表示重建距離z'= 168 μm、z'= 250 μm、z'= 382 μm 時RBC 重建圖像, 其中箭頭表示該重建距離下可清晰觀察到的單個RBC 樣品, 依次確定全息圖中185 個RBC 的深度位置, 聚焦RBC 中心位置的光強度具有局部最小值, 可通過光強度大小判斷每一個RBC 的橫向位置 (x,y) .

綜上所述, 可以獲取全息圖記錄的每一個RBC 的三維信息, 其相關x,y,z測量精度如表1所示, 由表可知RBC 的橫向信息與軸向信息的測量精度均較高, 能準確獲取RBC 的三維信息. 進一步地, 通過CMOS 像感器的連續拍攝, 記錄RBC 的運動軌跡全息圖, 通過同軸全息重建算法對RBC 進行重建, 分析FEP 微流管內RBC 的三維運動, 實現RBC 動態追蹤, 得到RBC 的實時空間分布.

圖6 RBC 的全息圖與重建圖像 (a) CMOS 像感器拍攝FEP 微管內RBC 所得全息圖; (b), (c), (d)分別表示不同重建深度下RBC 重建圖像, 箭頭表示聚焦的RBC[51]Fig. 6. The hologram and the reconstruction images of RBCs. (a) The RBC hologram obtained from CMOS; (b),(c), and (d) represent RBC reconstruction images at different reconstruction depths, respectively. Each focused RBC is shown by an arrow[51].

表1 RBC 橫向信息與軸向信息的測量精度[51]Table 1. The lateral and axial measurement accuracy of RBC[51].

3.3 離軸數字全息顯微系統識別RBC 形狀

為實現高分辨率和無損地識別不同形狀生物細胞[60,61], 韓國朝鮮大學的Moon 等[54]提出一套離軸數字全息顯微系統用于測量RBC 的三維體積, 該系統引入分水嶺分割算法, 它是一種基于拓撲理論的數學形態學分割算法, 將圖像看作一幅地形圖, 圖像中每一點像素灰度值表示該點海拔高度, 將每一個局部極小值及其影響區域作為集水盆, 集水盆的邊界形成分水嶺, 根據圖像自身灰度差異進行區域劃分[62]. 通過分水嶺分割算法對重建RBC 相位圖像進行區分, 單獨提取每一個RBC, 根據單個RBC 自身表面積和相位值使用自動算法分別計算每個RBC 的體積, 生成全部RBC 的三維體積分布. 為驗證該方案的可行性, 通過上述方案分別計算生成口腔狀RBC 與盤狀RBC 的三維體積分布, 實驗表明, 所提出的離軸數字全息顯微成像方法可以對各種形狀RBC 進行三維體積定量分析, 其光學系統的分析與實現可根據圖7 所述結構進行描述, 激光器發射682 nm 波長紅光, 由分束器(BS1)將入射光分為兩束, 第一束經過待測樣本發生衍射, 經顯微物鏡MO(40×/0.75NA)放大后作為放大物光先后經過反射鏡與分束器后入射至CCD 像感器, 第二束光作為參考光, 通過分束器(BS2)后與放大物光在CCD 像感器上發生干涉, CCD 像感器記錄生成的全息圖.

圖7 離軸數字全息顯微成像系統測量不同形狀RBC 三維體積[54]Fig. 7. Off-axis digital holographic microscopy system for measuring RBCs’ three-dimensional volume of different shapes[54].

根據離軸全息重建算法對記錄的全息圖進行重建生成RBC 相位圖, 單個RBC 的體積計算公式如下:

其中,N表示RBC 的像素總數,p表示像素尺寸大小,φ表示RBC 內每個像素對應的相位值,M為顯微物鏡的放大倍數. RBC 折射率nrbc=1.396 ,HEPA 介質折射率nm=1.3334 .

為了研究不同形狀RBC 的三維體積, 在圖7 中選取了口腔狀(stomatcyte)RBC 與盤狀(discocyte)RBC 進行對比. 口腔狀RBC 與盤狀RBC 重建相位圖像如圖8(a)、8(b)所示, 通過分水嶺分割算法去除相位圖像背景與RBC 的互相重疊部分, 對RBC 進行逐一區分, 如圖8(c), 8(d)所示, 可以看出口腔狀RBC 的重建相位圖像與盤狀RBC 的重建相位圖像中相鄰細胞間均明顯隔開, 能夠實現單個口腔狀RBC 或盤狀RBC 的提取. 為了對RBC的三維結構進行更合理地分析, 通過標記分水嶺分割算法對圖8(e)所表示的單個RBC 進行標記, 得到標記后的單個RBC 的整體C 部分, 如圖8(h)所示, 經過標記分水嶺分割算法計算得到單個RBC 的A 部分, 如圖8(f)所示. 通過減法運算, 將單個RBC 的整體C 部分減去A 部分可得到RBC的B 部分, 如圖8(g)所示.

圖8 RBC 的相位重建圖像 (a)重建的口腔形狀RBC 相位圖像; (b)重建的盤狀RBC 相位圖像; (c)重建后經分水嶺算法分割的口腔形狀RBC 相位圖像; (d)重建后經分水嶺算法分割的盤狀RBC 相位圖像; (e)分割的單個RBC 相位圖像; (f)(g)(h)分別經標記分水嶺算法進一步分割得到的單個RBC 的A、B、C 部分[54]Fig. 8. The reconstructed phase image for RBCs (a) The reconstructed phase image for RBCs having a stomatocyte shape; (b) the reconstructed RBCs phase image for RBCs having a discocyte shape; (c)the segmented phase image for RBCs having a stomatocyte shape; (d) the segmented phase image for RBCs having a discocyte shape; (e) the segmented phase image for single RBC(f), (g) and(h) represent the A, B and C parts by the marker-controlled watershed algorithm in RBC, respectively[54].

為了對兩種不同形狀RBC 的三維體積進行定量研究, 首先計算出最佳重建距離下RBC 的相位分布, 再根據(11)式分別對每個口腔狀或盤狀RBC 的A 與B 部分進行計算, 得到平均體積和均方根誤差的數值, 結果如表2 所示. 由表可知盤狀RBC 的A 部分體積大于口腔狀RBC 的A部分體積, 而B 部分體積小于口腔狀RBC 的B部分體積, 根據進一步的統計假設檢驗證明了RBC 的三維體積隨著其形狀的改變會出現較大區別[63]. 綜上所述, 根據RBC 的體積自動計算算法結合分水嶺分割算法可以研究不同形狀RBC 的三維體積并且實現實時檢測RBC 的三維體積變化, 對生物醫學領域中細胞成像發展研究具有極大價值.

表2 兩種不同形狀RBC 的A、B 部分的三維體積Table 2. The different shapes of RBC’s three-dimensional volume of A and B parts.

3.4 離軸數字全息顯微系統測量RBC 體積

數字全息的定量相襯成像在顯微生物成像中具備了多種新應用[64], 其中包括微觀對象的評估,如生物細胞生長觀測[65], 由滲透壓引起的細胞體積變化檢測[66], RBC 監測[67]. 這些研究中通常使用定性指標即通過重建圖像視覺質量來進行焦平面評估, 因此獲取結果具有不確定性[68?70]. 基于此, 美國杜克大學 Wax 等[53]設計了一套離軸數字全息系統并提出了一種數字重聚焦方法, 通過定量指標評估焦平面并描述微觀物體的方法, 在嘗試多種最佳焦平面的評估方法并對比每種方法的重建圖像效果后, 選擇振幅方差最小化作為定量指標,由于其在評估物體透射的相位圖像時, 在深度方向聚焦清晰度與聚焦范圍均具有較大優勢, 且計算量較小, 可找到理想的焦平面位置. 其光學系統的分析與實現可根據圖9 所述結構進行描述, 低相干激光器發射589 nm 黃光, 分束器將黃光分為兩束,由后向反射鏡實現兩束光的光程匹配, 第一束光通過后向反射鏡調整光程入射至RBC 發生衍射,RBC 置于樣品臺上, 樣品臺可移動以控制樣品與焦平面距離, 衍射光經過顯微物鏡MO( 40×/0.75NA)放大后作為放大物光入射至CMOS 像感器; 另一束光作為參考光通過后向反射鏡調整光程經顯微物鏡MO(40×/0.75NA)放大后作為放大參考光與放大物光在CMOS 像感器上發生相干疊加, 形成全息圖, 由CMOS 記錄生成的全息圖.

圖9 離軸數字全息顯微系統研究散焦現象對RBC 三維體積測量影響[53]Fig. 9. Off-axis digital holographic microscopy system for investigating the effect of defocus on RBC three-dimensional volume measurement[53].

根據離軸全息重建算法對記錄全息圖進行重建, 獲得RBC 重建圖像. 為檢驗該定量指標評估焦平面后重建物體的指標, 定義了微觀物體的光學體積, 公式如下

其中, 光程長度的變化 ?OPL(x,y)=?φ(x,y)/k=?n(x,y)·h(x,y)用于描述微觀物體,h(x,y) 表示物體的厚度(高度)圖, ?n(x,y) 表示物體折射率圖,k=2π/λ表示波數, 合并 ?OPL(x,y) 參數產生物體光學體積OV的度量, 若物體的折射率均勻,即 ?n(x,y)=?, 說明物體光學體積與物體真實體積線性相關,OV=?.

為進一步比較定性指標與定量指標分析中分別評估得到的焦平面效果, 根據焦平面各自對應的重建物體圖像的光學體積對比原物體光學體積真值, 以誤差值衡量所得焦平面的效果. 引入直徑為8 μm 的聚苯乙烯微球(直徑變異系數1%)作為樣品, 首先根據人工聚焦方法即重建圖像視覺質量來進行焦平面評估, 移動樣品臺通過CMOS 像感器記錄微球全息圖, 在不同重建距離下對微球圖像進行重建, 通過人工觀察重建微球的圖像質量確定最佳焦平面位置, 再計算重建微球的相位分布, 設定此時重建微球球心所在平面位置為z= 0. 根據提出的數字重聚焦方法計算z= ± 50 μm 范圍內微球球心處的振幅方差, 確定振幅方差最小值位置為z= –2.46 μm, 將其作為最佳焦平面位置, 并計算微球的相位分布. 通過(12)式分別計算人工聚焦方法與數字重聚焦方法所得微球光學體積, 對比微球光學體積真值(1%的直徑變異系數對應于3%光學體積變異系數), 結果如表3 所示. 表中±7.98 fL 表示聚苯乙烯微球3%的光學體積變異系數, 人工聚焦方法測得微球OV= 281.7 fL, 以5.90%比率高于微球真值OV= 266.0 fL, 而數字重聚焦方法測得微球光學體積為OV= 266.6 fL,僅以0.22%比例高于微球光學體積真值, 說明根據數字重聚焦方法所得焦平面更具可靠性, 將其認定為理論最佳焦平面, 此時人工聚焦方法所得z=0 處距離理論最佳焦平面距離為2.46 μm, 側面反映了即使在2.46 μm 的輕微散焦距離下, 所測得光學體積結果仍會產生明顯誤差. 為研究散焦距離變化對微球光學體積測量結果的影響, 通過移動控制臺改變微球散焦距離, 實時測量對應微球的光學體積, 對比最佳焦平面處測得微球的光學體積, 可以得到, 即散焦距離每改變1 μm,其光學體積測量結果相差2.2%, 當散焦距離達到± 10 μm 時, 會出現相位包裹偽影進一步降低測量精度, 而根據數字重聚焦方法可以解決散焦的問題, 進一步證明了該方法在待測樣本光學體積測量時的高精度優勢.

表3 人工聚焦方法與數字重聚焦方法測得微球光學體積對比[53]Table 3. Comparison of OV measured by manually-focused and digitally-refocused methods[53].

為顯示生物樣品中光學體積測量中散焦的影響, 選取單一RBC 作為測試樣品, 在約200 μm 范圍人工選取7 個不同焦平面, 分別對RBC 進行重建得到相對應的振幅圖與相位圖, 如圖10(a)中首行與中間行(A—G)所示.

圖10 對單一RBC 圖像進行數字重聚焦與相應的RBC 光學體積測量 (a)通過人工聚焦方法與數字重聚焦方法對單一RBC 重建所得振幅圖與相位圖; (b)A—G 的振幅方差分布; (c)RBC 在人工聚焦方法所得光學體積(黑線)與數字重聚焦方法所得光學體積(藍線), 光學體積OV 表示為平均值 ± 標準差[53]Fig. 10. Digital refocusing of a single red blood cell image and corresponding optical volume measurements. (a) The amplitude and phase images by the manually-focused method and digitally-refocused method from a single RBC; (b) amplitude variance metric of holograms A-G; (c) computed OV of RBC from manually-focused phase images(black) and digitally-refocused phase images(blue).OV reported as mean ± standard deviation[53].

其中, D 圖表示重建微球視覺質量最佳時所得微球的振幅圖與相位圖, D 圖的焦平面位置即為定性指標中最佳焦平面位置, 分別計算A—G 對應的振幅方差最小值可得到7 個焦平面相對定量指標中理論最佳焦平面位置(z=0 )各自的散焦距離,如圖10(b)所示, 人工聚焦方法所得重建RBC 圖像A—G 的光學體積計算結果如圖10(c)中黑色線所示, 隨著散焦距離的變化, RBC 在OV=5.15—7.03 fL (σ=13.55% )范圍內波動, 散焦距離較大影響了RBC 體積測量數值的可靠性; 而通過數字重聚焦方法獲取振幅方差最小值以得到理論最佳焦平面位置, 重建所得RBC 相位圖如圖10(a)末行所示, 通過數字重聚焦方法所得A—G 的相位圖相似, RBC 的光學體積僅在OV= 5.54—5.70 fL(σ=0.99% )范圍內波動, 進一步證明數字重聚焦方法對光學體積測量的實用性與可靠性. 綜上所述, 根據數字全息顯微成像方法表征微觀物體三維體積時, 重建圖像的散焦現象是潛在的重要誤差來源, 提出的數字重聚焦方法能夠以定量指標評估焦平面, 對重建RBC 進行光學體積測量, 實時獲得準確RBC 的三維體積, 可廣泛應用于人體的心血管疾病與糖尿病的病理研究之中.

3.5 光鑷輔助離軸數字全息系統

帕金森氏疾病(PD)作為一種中腦區神經元丟失的神經退行性疾病, 隨著PD 發病機理研究的不斷深入, 在各種病理假設之中, 氧化應激是引起細胞功能障礙與細胞凋亡的主要因素[43,44]. 在隨后的研究中, 證實了PD 中氧化應激的存在, 根據以往研究發現氧化應激還會誘導細胞變形[71], 因此,研究氧化應激下細胞的機械性能對PD 的研究具有重要意義. 基于此, 祝連慶等[72]設計了一套光鑷輔助離軸數字全息顯微系統, 用于研究氧化應激下RBC 的高度與三維體積變化, 有助于帕金森氏疾病的深入研究. 其光學系統的分析與實現可通過圖11 所述結構進行描述, 引入含有鏈霉親和素涂層的聚苯乙烯微球(直徑4 μm)以2∶1 的濃度比與人類RBC 樣本進行混合至磷酸鹽緩沖溶液之中(PBS, pH = 7.4), 篩選出RBC 水平兩側分別連結一個微球的樣品作為實驗待測對象, 如圖11中待測樣本放大圖所示. 具有統一參數的水平兩側微球通過光鑷控制一側微球靜止, 另一側恒定慢速移動以實現RBC 的拉伸, 考慮到RBC 的尺寸并非絕對均勻, 因此, 從滿足實驗條件細胞中隨機選取10 個RBC 取計算結果平均值, 通過如圖11 所示光學系統進行全息圖的記錄, 光學系統由離軸數字全息顯微系統與光鑷組成, 632.8 nm激光器發射紅光通過分束器(BS1)分為物光與參考光, 物光通過油浸物鏡MO(60× /1.40NA)放大后與參考光于CMOS 像感器記錄面上發生相干疊加, 形成干涉全息圖, 參考光的反射角度通過反射鏡的傾斜度調整, 用于控制干涉條紋的方向與寬度, 1064 nm激光器(5 W, Nd: YAG)作為光鑷的照明光源, 發射紅外光通過二向色鏡(DM)進入倒置顯微鏡的后孔形成光鑷, 其中聲光偏轉器(AOD)通過聲光的相互作用, 改變激光角度從而控制RBC 水平兩側微球的位置以實現拉伸.

利用離軸全息重建算法對CMOS 像感器所記錄的全息圖進行重建, 并通過重建相位圖像的參數計算不同氧化應激濃度下RBC 的高度與三維體積, 公式如下:

圖11 光鑷輔助離軸數字全息顯微系統測量RBC 體積[72]Fig. 11. Off-axis digital holographic microscopy system with optical tweezer for measuring RBCs’ three-dimensional volume[72].

其中,λ代表全息顯微系統中的激光波長(632.8 nm),φ代表相位值, RBC 的折射率ncell=1.37 ,溶液的折射率nmedium=1.34,V表示RBC 的體積,h(xi,yj)表示RBC 上任意點的 (xi,yj)的高度,M和N表示RBC 的橫向坐標范圍.將過氧化氫(H2O2)溶液(0.1 mol/L)在PBS中稀釋至50、100、150 和200 μmol/L 對照組, 置于5 個培養皿中, 與符合實驗條件附著有微球的RBC 樣品進行混合. 為具體分析陷阱的拉伸力與細胞變形之間的關系, 需對光阱剛度值進行校準,校準分為已知力校準方法和熱漲落方法[73,74]. 熱漲落是系統在平衡狀態下相對于其平均狀態的隨機偏差, 是系統溫度的基本體現, 根據Boltzmann 統計分析熱噪聲下隨機運動粒子的布朗軌跡, 可以獲得粒子位移概率分布圖, 接著通過計算粒子的電勢獲得光阱剛度k, 該值被校準為14.23 ± 0.46 pN/μm.校準后保持激光功率不變以確保光阱強度均勻. 光阱拉伸力F由光阱剛度k和微球與陷阱之間中心差 ?x共同確定, 控制陷阱拉伸力在0~3 pN 之間變化. 將RBC 所能達到的最大高度值定義為H,同時為了清楚觀察H的變化, 在圖12 中分別顯示了H大于2 μm 的高度部分, 隨著氧化應激濃度的增大, RBC 的細胞膜出現粗糙度且其邊緣呈現鋸齒狀. 分別計算五項氧化應激濃度下測量所得10 個RBC 的最大高度與三維體積的平均值, 能夠更直觀比較RBC 在不同濃度氧化應激下高度與三維體積的變化, 繪制的相關折線圖如圖13 所示,H/H0表示拉伸過程中細胞的最大高度H與未拉伸前RBC 的原始高度H0的比值, 圖13(a)表示0—200 μmol/L 氧化應激濃度下RBC 高度比H/H0與陷阱拉伸力的關系, 反映了RBC 的微形變機械性能, 由圖可知, 隨著氧化應激濃度的增加,RBC 的變形能力下降, 因此氧化應激下的RBC 其機械性能會受到嚴重影響. 圖13(b)表示0—200 μmol/L 氧化應激濃度下RBC 三維體積的變化, 由圖可知, 對比不加入H2O2的常規PBS 中的RBC, 氧化應激使RBC 的體積減小, 在50—200 μmol/L 的H2O2濃度下, RBC 的體積變為初始體積的90%—87.3%, 該結果與過往研究相一致[75]. 綜上所述, 在氧化應激下RBC 的機械性能會發生變化, 微形變能力下降, 而氧化應激是導致PD 的主要因素, 該光鑷輔助離軸數字全息顯微成像的技術將助于PD 的基礎研究與臨床應用.

圖12 0?3 pN 陷阱拉伸力變化下不同濃度(0?200 μmol/L)氧化應激下重建RBC 的高度變化, 顏色深淺代表RBC 高度的大小[72]Fig. 12. Height change of reconstructed RBCs under different concentrations of oxidative stress (0?200 μmol/L). Four images in each group are corresponding to trap force varying from 0?3 pN. Color bar represents different thickness[72].

圖13 RBC 在不同濃度氧化應激下的性能 (a)不同濃度氧化應激下RBC 最大高度H 與陷阱拉伸力關系; (b)不同濃度氧化應激下RBC 體積[72]Fig. 13. Performance of RBC under different oxidative stress. (a) The relationship between the maximum height H of RBC and the trap tensile force under different oxidative stress; (b) the volume of RBC under different oxidative stress[72].

4 結論與展望

本文簡要介紹了數字全息顯微成像由全息圖記錄到全息圖重建的原理, 分析討論了不同數字全息顯微系統在RBC 成像的應用, 在數字全息顯微系統方面, 主要分為同軸、離軸、光鑷輔助離軸的數字顯微全息成像系統. 同軸數字全息顯微成像系統光路簡單, 可采用激光作為光源, 但其對光源的相干性以及系統的穩定性要求相對離軸數字全息顯微成像系統較低, 可通過弱相干光源(如LED)代替激光作為光源, 采用弱相干光源避免了激光產生的散斑噪聲, 可以得到較好的生物細胞圖像, 降低成本, 裝置結構輕便, 便于集成. 離軸數字全息顯微成像系統引入參考光光路, 通過改變參考光入射角度, 實現重建像與孿生像的分離. 光鑷輔助離軸的數字全息系統中, 光鑷作為輔助檢測的手段,可給予活體細胞不同拉伸力使之產生形變, 并根據離軸數字全息系統對細胞進行實時成像, 進而對活體細胞的性能進行更深一步的研究. 在不同的實驗要求下, 可根據具體需求選擇合適的全息顯微成像系統以實現實驗目的.

在細胞三維信息提取的過程中, 數字重聚焦方法選擇圖像在不同距離下振幅方差最小值作為定量指標代替圖像視覺質量判斷的定性指標, 以評估RBC 的焦平面位置, 減少由于散焦造成的三維體積測量誤差. 或根據圖像清晰度量化算法也可以實現RBC 焦平面評估, 確定重建距離, 通過重建圖像信息進一步獲取RBC 表面積、相位值, 并根據入射光波長、折射率等參數信息計算RBC 高度與三維體積, 也可通過反向重建放大物光復振幅的實部信息來判斷RBC 幾何形狀的變化; 對于多個RBC 的情況, 引入分水嶺分割算法對重建RBC 相位圖像進行分割, 單獨提取每一個RBC 后, 再通過上述算法分別計算每個RBC 的高度與三維體積并取對應平均值獲得RBC 整體三維信息. 更進一步地, 通過光鑷控制RBC 拉伸產生微形變, 根據細胞的三維信息探究RBC 在不同濃度氧化應激下的形變能力.

通過數字全息顯微成像技術有效實現細胞三維體積的測量, 并且可以檢測出細胞幾何形狀的微小形變, 對比于當前顯微成像領域廣泛運用的激光掃描共聚焦顯微技術, 近場光學顯微技術, 光學相干層析成像技術, 其非接觸、無損性、定量化、實時性的優勢更有利于細胞生物學參數與細胞性能的研究. 未來數字全息顯微成像技術有兩個大方向可以進行持續探索, 一方面可進一步發揮先進光學器件優勢, 在硬件上實現分辨率的提升, 并設計更優化的重建算法, 通過更少的數據量與更精確的焦平面評估來實現更高分辨率, 更高精度的細胞成像;另一方面由于系統的分辨率受限于顯微物鏡數值孔徑, 可采用合理的數字全息無透鏡成像方法, 待測細胞與像感器之間不存在光學透鏡, 因此可避免透鏡表面灰塵造成的雜散光影響, 且簡化了成像的光路, 具有更大的成像視野, 同時滿足細胞顯微成像的大視場與高分辨率需求, 便于細胞成像研究方面的集成化運用. 同時, 基于深度學習的圖像重建算法的發展, 可根據增加像素數的方式在傳統方法的基礎上進一步提升細胞成像的分辨率, 提高成像質量.

主站蜘蛛池模板: 国产丰满成熟女性性满足视频| 久久综合激情网| 99re视频在线| 欧美特黄一级大黄录像| 久久国产精品波多野结衣| 99手机在线视频| 91在线高清视频| 国产SUV精品一区二区6| 呦视频在线一区二区三区| 激情国产精品一区| 3D动漫精品啪啪一区二区下载| 日本国产精品| 久草热视频在线| 亚洲男人的天堂久久香蕉网| 精品成人一区二区三区电影| 国产精品人莉莉成在线播放| 天天综合网站| 91尤物国产尤物福利在线| 国产精品成| 国产免费羞羞视频| 国产乱子伦一区二区=| 国产午夜无码片在线观看网站| 久久精品一品道久久精品| 久久婷婷综合色一区二区| 亚洲福利片无码最新在线播放| 欧美在线综合视频| 日韩精品一区二区三区视频免费看| 日韩区欧美区| 国产无码性爱一区二区三区| 欧美一级在线| www.日韩三级| 91久久青青草原精品国产| 国内精品91| 国内精自视频品线一二区| 影音先锋丝袜制服| 亚洲 欧美 日韩综合一区| 欧美日韩亚洲国产主播第一区| 毛片免费在线视频| 亚洲欧美色中文字幕| 国产黑丝一区| 亚洲天堂视频在线观看免费| 欧美激情成人网| 国产在线观看高清不卡| 天堂va亚洲va欧美va国产| 国产精品林美惠子在线播放| 亚洲视频免| 国产午夜精品一区二区三区软件| 精品成人一区二区三区电影 | 免费国产小视频在线观看| 熟女视频91| 看av免费毛片手机播放| 四虎精品国产AV二区| 国产一区二区人大臿蕉香蕉| 99热这里只有精品2| 2019国产在线| 影音先锋亚洲无码| 毛片免费高清免费| 国产成人无码AV在线播放动漫| 国产主播福利在线观看| 最新痴汉在线无码AV| 国产一级毛片网站| 97精品伊人久久大香线蕉| 国产精品男人的天堂| 国产精品伦视频观看免费| 不卡无码网| 亚洲A∨无码精品午夜在线观看| 毛片在线播放网址| 男女男免费视频网站国产| 91欧洲国产日韩在线人成| 国产精品思思热在线| 中文字幕无码制服中字| 色天堂无毒不卡| 国内a级毛片| 91在线无码精品秘九色APP| 亚洲美女久久| 午夜视频在线观看区二区| 亚洲精品制服丝袜二区| 久久国产精品影院| 亚洲男人的天堂网| 亚洲天堂网在线观看视频| 国产欧美日韩免费| 国产欧美日韩va另类在线播放|